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ABSTRACT 

A major goal in the control of complex mechanical systems such as spacecraft 
rocket engines, advanced aircraft, and power plants is to achieve high performance with 
increased reliability, availability, component durability, and maintainability. The 
current practice of decision and control systems synthesis focuses on improving 
performance and diagnostic capabilities under constraints that often do not adequately 
represent the materials degradation. In view of the high performance requirements of 
the system and availability of improved materials, the lack of appropriate knowledge 
about the properties of these materials will lead to either less than achievable 
performance due to overly conservative design, or over-straining of the structure 
leading to unexpected failures and drastic reduction of the service life. The key idea in 
this report is that a significant improvement in service life could be achieved by a small 
reduction in the system dynamic performance. The major task is to characterize the 
damage generation process, and then utilize this information in a mathematical form to 
synthesize a control law that would meet the system requirements and simultaneously 
satisfy the constraints that are imposed by the material and structural properties of the 
critical components. 

The concept of damage mitigation is introduced for control of mechanical 
systems to achieve high performance with a prolonged life span. A model of fatigue 
damage dynamics is formulated in the continuous-time setting, instead of a cycle-based 
representation, for direct application to control systems synthesis. An optimal control 
policy is then formulated via nonlinear programming under specified constraints of the 
damage rate and accumulated damage. The results of simulation experiments for the 
transient upthrust of a bipropellant rocket engine are presented to demonstrate efficacy 
of the damage-mitigating control concept 
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CHAPTER 1 


INTRODUCTION 

1.1. Motivation and Statement of the Research Problem 

A major goal in the control of complex mechanical systems such as spacecraft 
rocket engines, advanced aircraft, and power plants is to achieve high performance with 
increased reliability, availability, component durability, and maintainability [Lorenzo 
and Merrill (1991a); Ray et al. (1993a, 1993b)] which manifest the following 
requirements: 

• Extension of the service life of the controlled process; 

• Increase of the mean time between major maintenance actions; 

• Reduction of risk in the integrated control-structure-materials system design. 

Therefore, the control systems need to be synthesized by taking performance, mission 
objectives, service life, and maintenance and operational costs into consideration. The 
current practice of decision and control systems synthesis focuses on improving 
dynamic performance and diagnostic capabilities under constraints that often do not 
adequately represent the materials degradation. The reason is that the traditional design 
is based upon the assumption of conventional materials with invariant characteristics. 
In view of high performance requirements and availability of improved materials that 
may have significantly different damage characteristics relative to conventional 
materials, the lack of appropriate knowledge about the properties of these materials will 
lead to either of the following: 

• Less than achievable performance due to overly conservative design; or 

• Over-straining of the structure leading to unexpected failures and drastic reduction 
of the useful life span. 

For example, reusable rocket engines present a significantly different problem in 
contrast to expendable propulsion systems that are designed on the basis of 
minimization of weight and acquisition cost under the constraint of specified system 
reliability. In the reusable rocket engines, multiple start-stop cycles cause large thermal 
strains; steady-state stresses generate inelastic strains; and dynamic loads induce high 
cyclic strains leading to fatigue failures. The original design goal of the Space Shuttle 
Main Engine (SSME) was specified for 55 flights before any major maintenance, but 
the current practice is to disassemble the engine after each flight for maintenance 
[Lorenzo and Merrill (1991a)]. Another example is design modification of the F-18 
aircraft as a result of conversion of the flight control system from analog to digital, 
which would lead to a significant change in the load spectrum on the airframe structure. 
In this case, a major goal of the vehicle control systems redesign should be to achieve a 
trade-off between flight maneuverability, and durability of the critical components 
[NoU etaL (1991)]. 

As the science and technology of materials continue to evolve, the design 
methodologies for damage-mitigating control must have the capability of easily 
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augmentation of current system-theoretic and Al-based techniques for synthesis of 
decision and control laws with governing equations and inequality constraints that 
would model the mechanical properties of the materials for the purpose of damage 
representation and failure prognosis. The major challenge in this research is to 
characterize the damage generation process, and then utilize this information in a 
mathematical form for synthesizing the control algorithms in complex mechanical 
systems. 

Although a significant amount of research has been conducted in each of the 
individual areas of Control Engineering and Analysis & Prediction of Materials 
Damage , integration of these two disciplines has not received much attention. 
Recently, Noll et al. (1991) have pointed out the need for interdisciplinary research in 
the fields of active control technology and structural integrity, specifically fatigue life 
analysis, in view of the integrated structural and flight control of advanced aircraft. 
Lorenzo and Merrill (1991b) have proposed a concept of damage mitigation and failure 
prognosis in the framework of reusable rocket engines for space propulsion in the 
context of intelligent control and diagnostics. This intelligent control system is 
hierarchically structured in which the top level coordinates the major functionalities of 
life extending control, adaptive control, real-time diagnostics and prognostics, and 
sensor/actuator fault accommodation. A major goal is to ensure safety and achieve the 
mission objectives while arbitrating the potentially conflicting requirements of 
optimum performance and structural durability. However, this concept of damage 
mitigation is not restricted to intelligent control and diagnostics, and can be applied to 
any system where structural durability is an important issue. 

1.2. Contributions of the Report 

The research reported in this report addresses the two different disciplines of 
Control Systems Engineering and Fracture Mechanics. The major contribution of this 
research is the development of a usable damage model and its application to the damage 
mitigating control of complex mechanical structures. In contrast to the usual notion of 
expressing the fatigue damage rate relative to the number of cycles, a concept of time 
derivative of the damage has been developed based on the conventional life prediction 
methods. A unique advantage of this damage model in the continuous-time setting is 
that it can be directly incorporated within the control system structure to provide the 
necessary information for on-line damage-mitigating control as well as for off-line 
synthesis of the control system. 

The above concept of continuous-time damage modeling can be related to a 
variety of cycle-based life prediction methods. This is essential for synthesis of a 
damage-mitigating control system since engineering applications may adopt different 
methodologies for damage analysis due to the diversity of engineering materials and 
damage accumulation mechanisms. 

1.3. Organization of the Report 

The report is organized in six chapters, including the introduction, and five 
appendices. Chapter 2 presents a general structure of the damage-mitigating control 
system along with a description of each component Chapter 3 describes the procedure 
of synthesizing an open-loop control law via nonlinear programming to optimize the 
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system performance without violating the imposed constraints. Chapter 4 reviews the 
current life prediction methods and derives die continuous-time damage model based 
on an experimentally verified approach. The control system is then simulated on the 
basis of a simplified model of the space shutde main engine (SSME), and the results of 
simulation experiments are discussed in Chapter 5. Chapter 6 summarizes the potential 
problems and difficulties in the current approach and suggests possible solutions and 
the areas of future research. Five appendices provide die necessary information to 
complement the results reported in the main body of this report Appendix A depicts 
the structure of a closed-loop control system with discussions on robustness and 
reliability. Appendix B explains the basic principle of a cycle counting approach for 
fatigue life prediction under spectral loading. Appendix C presents a mathematical 
model of a bipropellant rocket engine which is the plant under control. Appendix D 
describes the stress analysis of a turbine blade which is considered to be a critical 
component of the plant Finally, an approach for modeling of integrated fatigue- 
corrosion-creep damage is outlined in Appendix E. 



CHAPTER 2 

DAMAGE-MITIGATING CONTROL SYSTEM 


The damage-mitigating control system (DCS), also referred to as Life Extending 
Control system [Lorenzo and Merrill (1991b)], is intended to function independently or 
as an integral part of a hierarchically structured control system. The DCS may be 
centralized or distributed depending on the spatial location of the critical plant 
components and the overall objective of the mission. As stated in Chapter 1, the major 
challenge in the DCS design is to characterize the damage generation process such that 
this information can be directly applied to synthesize the control law. Therefore, this 
report focuses on the fundamental issue of formulating a generic structure of the DCS 
with the objective of achieving high performance while simultaneously maintaining the 
accumulated damage and the damage rate of the critical plant component(s) within the 
prescribed limits. 

2.1. General Structure of the Damage-Mitigating Control System 

Figure 2.1 shows a conceptual view of the DCS which consists of the following: 
(i) a physical plant under control; (ii) a damage predictor which models the structural 
and damage dynamics of the critical plant components; and (iii) a feedback controller 
which coordinates the plant states and damage mitigation. Referring to the control 
systems structure in Figure 2.1, the DCS synthesis consists of the following major 
tasks: 

• Task 1 : Formulation of the plant model to represent the dynamic characteristics of 

the physical process; 

• Task 2 : Formulation of the damage predictor to represent the dynamic characteristics 

of the structural and material properties of the critical plant components ; 

• Task 3 : Synthesis of an open-loop control policy via optimization of a cost 

functional (representing the mission objectives) under the constraints 
imposed by the plant model and the damage prediction model; and 

• Task 4 : Synthesis of a closed-loop control policy via feedback of the plant states and 

damage information to steer the plant along the open-loop trajectory without 
exceeding the prescribed damage constraints. 

As shown in Figure 2.1, the sequence of open-loop control commands, serve as 
the feedforward input to the plant to fulfill the mission objectives under the specified 
damage constraints. The plant model is a finite-dimensional state-space representation 
of the system dynamics (e.g., thermal-hydraulic dynamics of the space shuttle main 
engine or propulsion and aerodynamics of an aircraft). The estimated plant states and 
plant outputs are the inputs to the structural model which, in turn, generates the 
necessary information for the damage model. The output of the structural model is the 
load vector which may consist of (time-dependent) stress, strain, temperature, wear, 
level of corrosion in gaseous and aqueous environments, and other physico-chemical 
process variables at die critical point(s) of the structure. The damage model is a 
continuous-time as opposed to a cycle-based representation of life prediction so that it 
can be integrated with the plant and structural models for DCS synthesis. This damage 
model includes the effects of damage rate and accumulated damage at the critical 
point(s) of the structure that are subjected to the time-dependent load. The damage 
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state vector v(t) could indicate the level of micro-cracking, macroscopic crack length, 
wear, creep, corrosion, density of slip bands, etc., at one or more critical points, and its 
time derivative v(t) indicates how the instantaneous load is affecting the structural 



Figure 2.1 General Structure of Damage-Mitigating Control System 

components. These two vectors, v(t) and v(t), are the damage variables to be 
constrained in the synthesis of an optimal open-loop control policy, and also provide 
the important information for damage mitigation in the synthesis of the feedback 
control policy. The overall damage D(t) is a scalar measure of the combined damage at 
one or more critical points resulting from different effects (e.g., fatigue, creep, 
corrosion, or wear) relative to the preassigned allowable level v re f of the damage vector. 
Although D(t) may not directly enter the feedforward or feedback control loop, it can 
provide useful information for intelligent decision-making such as damage prognosis 
and risk analysis. In the closed-loop control, the state and damage feedback controller 
monitors and controls the plant and the damage states which may deviate from the 
desired values due to modeling uncertainties, external disturbances, and sensor noise. 
The role of this controller is to manage the possible conflicts between the plant 
performance and damage, and generate the appropriate feedback control efforts, u^, to 
maintain the plant response as close as possible to the nominal trajectory of the open- 
loop control policy. 
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22 . Formulation of the Damage-Mitigating Control System 

The dynamics of each block in Figure 2.1 need to be modeled in an appropriate 
form so that it can be direcdy applied to the synthesis of both the open-loop and closed- 
loop control laws. The plant dynamics and damage predictor are modeled by nonlinear 
(and possibly time-varying) differential equations which must satisfy the local Lipschitz 
condition [Vidyasagar (1992)] within the domain of the plant operating range. The 
structural model consists of solutions of dynamic equations representing the structural 
properties and the (mechanical and thermal) load conditions. These equations may be 
generated using a variety of analytical techniques ranging from finite element analysis 
to simple relationships of uniaxial stress and strain. However, the transformations 
relating the structural responses (such as nodal displacements) to the stress-strain 
behavior need to be included in the model. Since control systems are analyzed in the 
continuous-time setting, it would be necessary to express the damage equations in terms 
of time as the independent variable instead of the number of cycles. The continuous- 
time damage model is a necessity for practical implementation in the control synthesis 
procedure and, most importantly, generation of the intra-cycle damage information. 
The general formulation of the plant and damage dynamics and their constraints are 
represented in the deterministic setting as follows: 


Task period: Starting time to to final time tf 


Plant dynamics : 

dx 

— = f(x(t), u(t), t); x(to) = xo 
dt 

(2.1) 

Damage dynamics : 

^ = h(v(t), q(x,t), t); v(to) = vo; h> 0; Vt 
dt 

(2.2) 

Damage measure : 

D(t) = 5(v(t), v,ef) and D(t) e [0, 1] 

(2.3) 

Damage rate tolerance : 

0 £ h(v(t), q(x,t), t) < P(t); Vt s [to, tfl; 

(2.4) 

Accumulated damage tolerance : [v(tf) - v(to)] < <P 

(2.5) 


where 

x e R n is the plant state vector, 

u € R m is the control input vector, 

v e R r is the damage state vector; 

Vjref e K is the preassigned limit for the damage state vector; 

P(t) € R r and <p e R r are specified tolerances for the damage rate and accumulated 
damage, respectively; 

ye R p is the plant output vector; 

xe R n is an estimate of the plant state vector, 

q e R p is the load vector, and 

D e [0, 1] is a scalar measure of the accumulated damage. 
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The vector differential equations (2.1) and (2.2) become stochastic if the 
randomness of plant and material parameters is included in the models, or if the plant is 
excited by discrete events occurring at random instants of time [Sobczyk and Spencer 
(1992)]. The stochastic aspect of damage-mitigating control is a subject of future 
research. 

The proposed DCS in Figure 2.1 uses the concept of conventional state 
estimation and state feedback together with the information of damage rate and 
accumulated damage to generate an appropriate feedback control signal. Although this 
additional information renders the control system to be highly nonlinear, the dynamic 
performance and service life of the plant can be better managed with damage 
mitigation. Without the damage feedback, the controller depicted in Figure 2.1 would 
reduce to a conventional output feedback controller. This report focuses on the 
development of a continuous-time damage model and formulation of an open loop 
control policy. Although the synthesis of a closed loop control system is not within the 
scope of this report, this topic is briefly discussed in Appendix A. 



CHAPTER 3 

FORMULATION OF AN OPTIMAL OPEN-LOOP CONTROL POLICY 


Given an initial condition, an open-loop policy for the plant control is generated 
via nonlinear programming [Luenburger (1984)] by optimizing the plant dynamic 
performance without violating the prescribed upper bounds of the damage rate and the 
accumulated damage as discussed in Chapter 2. The open-loop control synthesis 
procedure is developed in this chapter. 


3.1. Problem Formulation 

The problem is to identify an open-loop control policy by minimizing a cost 
functional under the constraints of the damage rate and the accumulated damage over a 
period of time. The cost functional, J, is chosen to be the square of the weighted L2- 
norm of the plant states, damage rate, and control inputs. Dining the task period, [to, 
tf], the plant is steered from its initial state x(to) to the final state x(t£) along an optimal 
trajectory. In the formulation of an optimization problem, the control functions, u(t), t 
e [to, tf], needs to be discretized in time as a sequence [uk] so that the number of 
variables to be optimized is finite for implementation on a finite-state machine. 
Partitioning the task period [to, tf] into N steps at the discrete time instants {tk}, 
k=0,l,2,—, N, integration of the plant and damage dynamics in eqs. (2.1) and (2.2) of 
Chapter 2 yields the following results: 


Discretized plant dynamics : Xk+i = Xk + J t ^ +1 f(x(t),u(t),t) dt ; (3.1) 

Discretized damage dynamics : Vk+i = Vk + f t tk+1 h(v(t),q(x,t),t) dt; (3.2) 

l k 

where the plant state, damage state and damage rate vectors are denoted as: 

*k=x(tk); vk=v(tk); and v k =h(vk, q(xk, tk), tk) 

Then, the optimization problem is formulated in a general form as follows: 

N-l 

Minimize: J = X J k(*k«*k>“k) (3.3) 

k=0 

Subject to : 0 < h(vk, q(xk, tk), tk) < p(k) for k=l,2,3, •••, N; and 

(v N - v 0 ) < q> (3.4) 


where Xk = *k ” x ss and Ok = Uk “ u ss are deviations of the plant state vector and the 
control input vector from the respective final steady state values of x$s and u$s; and P(k) 

€ R r and 9 e R r are specified tolerances for the damage rate and accumulated damage, 
respectively. 
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The open-loop control law is synthesized by minimizing the cost functional in 
eq. (3.3) unden (i) the above inequality constraints in eq. (3.4); and (ii) the condition 
that, starting from the initial conditions x(to) and v(to), the state trajectory must satisfy 
the plant dynamic model in eq. (2.1). The design variables to be optimized are the 
discrete control inputs k=0,l,2, •••, N-l, and is not required to be identified 
because the optimization procedure is terminated at tN=tf. The flow chart shown in 
Figure 3.1 depicts a procedure for synthesizing an open-loop control policy via 
nonlinear programming. 

The first step in the nonlinear programming is to formulate a mathematical 
model of the plant dynamics in the state space form of eq. (2.1) with proper initial 
conditions. Based on the responses of the plant at the normal operating condition, the 
critical component(s) in which the damage is most likely to occur is (are) identified 
either by stress analysis or from the history of plant operations. Upon identification of 
the critical points for damage prediction, the structural modeling and the continuous- 
time damage modeling are used to formulate the damage dynamics in the form of eq. 
(2.2). The resulting damage rate and accumulated damage are to be constrained in the 
optimization procedure via nonlinear programming. The upper bounds of the 
constraints on damage rate and accumulated damage need to be appropriately selected 
by taking the mission objectives, the time interval between maintenance, service life 
and allowable risk into consideration. The initial damage is important due to its 
significant effects on the dynamics of nonlinear damage accumulation. The 
performance cost functional in eq. (3.3), which is minimized under nonlinear 
programming, is a function of the plant states, and damage rate vectors, representing a 
trade-off between system performance and damage. The weights are assigned to the 
plant states or selected output variables reflecting their relative impact on the system 
performance. If the damage constraints are appropriately chosen, it may not be 
necessary to include the damage rate in the cost functional. Upon selection of the cost 
functional and damage constraints, the task is to find an optimal control sequence {ut) 
in discrete steps from time to to tf. This optimal control sequence is then tested via 
system simulation and structural analysis to verify the plant performance and damage. 
If the results are satisfactory, the synthesis of an optimal policy for open-loop control is 
completed. Otherwise, the damage constraints should be revised and the optimization 
procedure is repeated to find a new control sequence (uk). It is also possible that an 
optimal solution may not exist due to overly stringent damage constraints. In that case, 
modifications of the constraints p or (p, or of the cost functional J are needed. If none 
of these revisions are advisable, revision of the damage model or the plant model 
should be considered. For example, in some situations, an alternative approach such as 
the fatigue crack growth model may be more suitable than the cyclic strain approach for 
damage modeling (see Section 4.6). 

In this report, a general purpose nonlinear programming package, namely the 
IMSL subroutine DNCONF [IMSL manual], has been adopted to solve the nonlinear 
optimization problem. Other nonlinear programming packages such as NPSOL [Gill et 
al. (1991)] of Stanford University are being considered for improving the computational 
efficiency in the future research. Following the optimization procedure in Figure 3.1, 
simulation experiments were carried out to evaluate the efficacy of the proposed 
damage-mitigating control. The results of simulation experiments are presented and 
discussed in Chapter S. 
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Figure 3.1 Nonlinear Programming Procedure 




11 


3 2 . Program Structure 

The software for open-loop control synthesis in simulation experiments is coded in 
Fortran using the IMSL subroutine DNCONF, and the program structure is shown in 
Figure 32 . The namelist data file on the left-hand top comer of Figure 3.2 contains the 
values of (i) the parameters of the plant model and the structural model, and (ii) 
material properties and damage constraints of the fatigue damage model. It serves as a 
data bank for the main program and supplies the necessary data to the subroutines 
through the COMMON statements. The MAIN program formats the output of the 
nonlinear programming (i.e., the optimal solutions: x ff , u ff , and v ff ), and these 
formatted data are to be used in the closed-loop control system in Figure 2.1 as the 
reference inputs. 


xff, uff, vff to closed-loop 



Figure 3.2 Program Structure for Open-Loop Control Policy via Nonlinear 
Programming 

The IMSL subroutine DNCONF, called by the MAIN program, uses a line 
search algorithm that minimizes the cost functional J. An user-supplied external 
subroutine FCN is called by DNCONF to calculate the objective functional and the 
values of constraints. The variables to be evaluated are the (open-loop) control input 
commands: u ff (0), uff(l), u ff (N-l), and each of these control commands is a vector 
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of dimension m. Therefore, the total number of the variables to be optimized is equal to 
mxN. At each step of the line search procedure, the sequence {u ff } is updated 
representing the improvement of the cost functional. The subroutine PLANT consisting 
of the plant dynamic model equations is called by FCN to calculate the state vector x ff 
for each input sequence, and the results are available to the structural model for 
evaluation of the load vector, q, which is the input to the damage model for prediction 
of the accumulated damage and damage rate. The damage rate, v ff , along with the state 
vector, x ff , and control input, u ff , are used by the subroutine COST to calculate the cost 
functional which, in turn, is fed back to die subroutine FCN for evaluation of the 
damage constraints. The line search continues until the cost functional converges and 
an optimal solution is reached. 



CHAPTER 4 


MODELING OF FATIGUE DAMAGE DYNAMICS 

Damage of mechanical structures is usually a result of fatigue, creep, corrosion, 
and their combinations [Suresh (1991)]. While the fatigue damage is cycle-dependent, 
the creep damage and corrosion damage are time-dependent. The prime focus in this 
report is on the representation of fatigue damage in a continuous-time setting. As 
discussed earlier, a time-dependent model of damage dynamics, having the structure of 
eq. (2.2), is necessary for analysis and synthesis of DCS. From this perspective, a 
dynamic model of fatigue damage has been formulated in the continuous-time setting. 
Although this damage model has a deterministic structure, it can be recast in the 
stochastic setting to include the effects of both unmodeled dynamics and parametric 
uncertainties. Augmentation of the fatigue damage model to include dynamics of 
corrosion and creep damage is briefly discussed in Appendix E. 

Because of the wide ranges in mechanical properties of materials, extensive 
varieties of experiments have been conducted for fatigue analysis, and many models 
have been proposed for fatigue life prediction in aerospace and ground vehicles 
[Bannantine et al. (1990); Dowling (1983)]. Each of these models expresses the 
damage dynamics by an equation with the number of cycles as the independent 
variable. In contrast, the damage dynamics in eq. (2.2) are expressed as a vector 
differential equation with respect to time, t, as the independent variable. The 
advantages of this representation are that it allows the damage model to be incorporated 
within the time-based structure of the constrained optimization problem and that the 
damage accumulated between any two instants of time can be derived even if the stress- 
strain hysteresis loop is not closed. The damage information at any desired point within 
a cycle is computed following the proposed approach. This concept is applicable to 
different models of damage dynamics such as those based on local strain or crack 
propagation. The strain-life and linear elastic fracture mechanics approaches for fatigue 
life prediction are briefly described below. 

• Cyclic Strain-Life : In this approach, the local stress-strain behavior is analyzed at 
certain critical points where failure is likely to occur. The local strain is directly 
measured from a strain gauge, or analytically computed. The local stress is 
estimated from the cyclic stress-strain curve. A cycle-based approach is then used 
to estimate the fatigue damage from the strain-life curves which are obtained 
through experiments at different levels of stress and strain amplitudes. The total 
accumulated damage is computed using the Palmgren-Miner rule [Miner (1945)] 
and subsequently modified via the damage curve approach [Bolotin (1989)]. 

• Linear Elastic Fracture Mechanics ( LEFM ) : The LEFM approach is built upon the 
concept of a physical measure of damage in terms of the crack length and the size of 
the plastic zone at the crack tip. The accumulated damage is computed by 
integrating the crack growth rate over the number of cycles. This is based on the 
crack growth rate equation being approximated by an exponential function of stress 
intensity factor range of the component [Bannantine et al. (1990)]. The component 
is assumed to fail when the crack reaches the critical length which, in turn, is 
determined from the fracture toughness of the component on the basis of 
experimental data. 
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Either of the above two approaches can be adopted for developing a continuous-time 
damage model depending on the availability of material data and the applications. The 
primary approach adopted in this report is based on cyclic strain-life. The next section 
reviews the governing equations and a procedure for fatigue damage prediction via the 
cyclic strain-life approach. An alternative approach based on LEFM is also presented 
in this chapter to demonstrate the applicability of the derivative concept 

4.1. A Review of Strain-Life Approach 

The fatigue damage is primarily controlled by the local strain at the critical 
point(s) of the component and is estimated by the strain-life curve which is essentially 
a plot of strain amplitudes (Ae/2) versus the number of cycles to failure (Nf) as 
obtained from constant amplitude fatigue tests on axially loaded specimens. A typical 
strain-life curve is shown in Figure 4.1. The relationship between Ae/2 and Nf can be 
expressed in the following form: 

-(2Nf) b + e f ' (2N f ) c (4.1) 

where 

Of : fatigue-strength coefficient 
b : fatigue-strength exponent 
£f : fatigue-ductility coefficient 
c : fatigue-ductility exponent 
E : elastic modulus 



Figure 4.1 Strain-Life Curve Showing Elastic and Plastic Components 
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The above coefficients are considered to be material constants which are available in 
the open literature [Tucker (1974); Boiler and Seeger (1987)] for steel and aluminum 
alloys. The constant amplitude loading experiments are normally conducted under 
completely reversed strain control with zero mean stress. However, different mean 
stress levels have been found to have important effects on the fatigue life. One method 

[Fuchs and Stephens (1980)] is to replace Of by o f ' - a m in eq. (4.1) to account for the 
mean stress effect such that 


Ae 

2 


Pf'-Pm 

E 


(2N f ) b + e f '(2N f ) c 


(4.2) 


where o m is the mean stress. Dowling (1983) further modified equation (4.2) as: 
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(4.3) 


The total number of cycles to failure (N f ) can then be computed by solving any of the 
above equations provided that the information of the stress and strain is known from the 
load history. The profiles of stress and strain, which are generated from either strain- 
control or stress-control experiments, form a hysteresis loop due to irrecoverable plastic 
strain under constant amplitude loading. The cyclic stress-strain curve is plotted by 
connecting the tips of hysteresis loops at different load amplitudes as shown in Figure 
4.2, which can be expressed in the following equation [Bannantine et al. (1990)]: 
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where n' : cyclic-strain hardening exponent which can be expressed as — ; 
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gf* . 


K' : cyclic-strength coefficient which can be expressed as 


-^7 is the elastic strain amplitude 
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is the plastic strain amplitude — — 
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(4.4) 
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Stress 



Figure 4.2 Cyclic Stress-Strain Curve 


For an unnotched component, the stress range A a can be computed from the 
dynamic system analysis. The strain range and total cycles to failure are then obtained 
from eqs.(4.1) to (4.4). In most notched components, however, a simple stress analysis 
at the notch is usually not possible. Instead, the local strain is first obtained by direct 
measurement at the notch if possible or by finite element analysis when it is difficult to 
measure due to complex geometry of the component. The strain is then converted to 
stress following the cyclic stress-strain characteristics of the specific material. 

The structural components under plant operations would usually be subjected to 
complex loading with varying amplitudes and frequencies. In that case, it becomes 
more difficult to recognize the cycles since the stress-strain hysteresis loops under 
spectrum loading are not closed on a cycle-by-cycle basis. To solve this problem, 
several approaches such as range-pair and rainfiow [Fuchs and Stephens (1980)] have 
been proposed to identify the cycles within the spectrum of the complex load history. 
All of these approaches attempt to extract the small cycles from the load history without 
losing the large cycles which significantly contribute to the damage. Since the rainfiow 
cycle counting is widely used in the strain-life approach, it has been adopted in this 
report for development of a continuous-time damage model. The concept of rainfiow 
cycle counting is briefly explained in Appendix B. Once a cycle is clearly identified, 
the linear damage of one cycle corresponding to a specific strain amplitude and mean 
stress is defined as: 



(4.5) 
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where N is the total number of cycles to failure at this stress and strain level. The 
accumulated damage of the load history is assumed to a linear summation of the 
individual cycles by the Palmgren-Miner's rule [Miner (1945)] as given below : 


n 1 

D=s ir 

i=l N fi 


(4.6) 


The fatigue damage is usually assumed to be independent of the frequency of 
the load and the shape of cyclic stress-strain hysteresis loop. However, when the 
component is exposed to high temperature, other effects such as creep may have to be 
taken into account and combined with the fatigue for damage computation. The 
modeling of fatigue-corrosion-creep is briefly discussed in Appendix E. 

4.2. Continuous-Time Damage Model 

This section introduces the development of a continuous-time damage model 
based on the strain-life approach as described above. The linear damage accumulation 
model in eqs. (4.5) and (4.6) is extended to detine the damage increment between any 
two points within a cycle as explained below: 

Following Figure 4.3, let the point O be the reference point of a reversal which 
is determined from rainflow cycle counting. Let A and B be two consecutive points on 
the same rising reversal, and let Na and Nb represent the total number of cycles to 
failure associated with constant amplitudes, OA/2 and OB/2, respectively. Denoting 
the damage of a rising reversal as S, the damage increment between the points A and B 
is defined as: 


AS = 


1 

n b 


1 

n a 


(4.7) 



Cycles to failure for constant 
amplitude OA : Na 


Cycles to failure for constant 
amplitude OB : N B 

Damage between A and B is defined by : 

1 1 
N b " N a 


Figure 4.3 Definition of Damage Between Two Points Within a Reversal 
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In eq. (4.7), it is assumed that the damage occurs only on the rising reversal, i.e., when 
the stress is monotonically increasing, and no damage is incurred during unloading, i.e., 
when the stress is monotonically decreasing. The assumptions are consistent with the 
physical phenomena observed in the fatigue crack propagation process [Suresh (1991)]. 
Given that A o is the stress increment between point A and point B, the average damage 
rate with respect to this stress change is equal to AS/Aa. Let At be the time interval 
from A to B, the average damage rate in terms of the stress can be transformed into the 
time domain by A8/At = (A8/Aa)x(Aa/At). Making At infinitesimally small, the 
instantaneous damage rate becomes 

® = Jim — x— = — x— (4.8) 

dt At— >0 Ac At da dt 

where the instantaneous stress rate da/dt can be generated either from direct 
measurements of strain rate or from finite element analysis, and dS/da is derived from 
the existing cycle-based formulae [Bannantine et al. (1990); Dowling (1983)]. 

In this report, Dowling's life prediction formula of eq. (4.3) and the stress-strain 
curve of eq. (4.4) are used to evaluate dS/da in eq. (4.8). Replacing Nf by 1/8, eq. (4.3) 
can be written as 


where e r is the total strain corresponding to the reference stress o r at the starting point 
of a given reversal as determined from the rainflow cycle counting method; and 

le ~ e r l/2 is the strain amplitude between the current point and the reference point The 

above equation does not provide a closed form solution for the predicted damage 8. 
The general approach to solve this problem is to separate eq. (4.9) into two different 
modes. The first term on the right side corresponds to the so called elastic damage 
mode 5e and the second term corresponds to the so called plastic damage mode Sn. 
These two damages, 8g and Sp, can be derived in an explicit form. For high cycle 
fatigue, the elastic damage usually yields more accurate prediction than the plastic 
damage, and vice versa for low cycle fatigue. Generally speaking, the transition 
between high cycle fatigue and low cycle fatigue is defined as point of intersection of 
the elastic strain-life curve and plastic strain-life curve as shown in Figure 4.1. In this 
report, it is proposed that the predicted damage 8 should be obtained as a weighted 
average of 8© and Sn where the weights depend on the relative accuracy of the elastic 
and plastic modes of damage computation under the current load condition. 

The elastic strain, £ re , and the plastic strain, Erp, corresponding to the reference 
stress. Or, are defined as: 



and Erp — Er “ Ere 


(4.10) 
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Eq.(4.9) can be split into the elastic damage mode and plastic damage mode as: 



(4.11a) 

(4.11b) 


where the elastic and plastic strain amplitudes, lEe-Erel/2 and l£p-£ml/2, respectively, are 
related to the state of stress following the cyclic stress-strain characteristics of the 
material as described in eq. (4.4) : 
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(4.12b) 


From eqs. (4.11) and (4.12), the closed form solutions for 5 e and 8p can be obtained in 
terms of stress instead of strain as given below: 
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(4.13a) 


(4.13b) 


Step changes in the reference stress o r can occur only at isolated points in the 
load spectrum. Since the damage increment is zero at any isolated point, the damage 
accumulation can be evaluated at all points excluding these isolated points which 
constitute a set of zero Lebesgue measure [Royden (1988)]. Exclusion of the points of 
step changes in o r does not cause any error in the computation of damage, and dOj/dt 
can be set to zero because o r is piecewise constant. Furthermore, since it is assumed 
that no damage occurs during unloading, the damage rate can be made equal to zero 
when a<o r - Then, the elastic damage rate dffe/dt and the plastic damage rate d5p/dt are 
computed by differentiating eqs. (4.13a) and (4.13b) as follows: 
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If 0^0*, then 


^ = 2-^ 


dt 


do 


n 


C-Cr 


2(°f -Om), 


do . 
— and 
dt 


d5.. 

2 

dt do 


'p d 


n 


1 

i |" o-o r y 
E f ' { 2 K J 




c 


'm 




do 

dt’ 


else dSe/dt=0 and d5p/dt =0. 


(4.14a) 


(4.14b) 


The damage rate d5/dt is obtained as the weighted average of the elastic and plastic 
damage rates such that 


d5 d5 e 

— = w — - 
dt dt 


d8 n 

+ (1 — w) — ~ 
v ' dt 


(4.15) 


where the weighting function, w, is selected on the basis of the elastic and plastic strain 
amplitudes in eqs. (4. 1 la) and (4.1 lb) as: 


w 


E e g tc 

e-e r 


and l-w= Ep £fp 

E-Er 


(4.16) 


Eqs. (4.14) to (4.16) are then used to obtain the damage rate at any instant The 
damage increment between two consecutive points tt and tk+i on the same reversal can 
be calculated by integrating the damage differential d8. The procedure for evaluation 
of the accumulated damage is described below: 


Assuming that there is no reference point change between tk and tk+i, i.e., no 
small cycle is closed between tk and tk+i, the damage increment during this interval is 
given as: 




f(l-w) 


ay 

dt J 


dt 


- w ( 8e(o(tk+l), Or) - 8e(o(tk), Or) ) + (1-w) ( 5p(0(tk+l), Or) - 5p(0(tk), Or) ) (4.17) 


where the weighting function is evaluated at the final point tk+i such that w = (Ee(tk+l) _ 

Ere)/(£(tk+l)“E r ). Furthermore, w is assumed to remain constant in the time interval (tk, 
tk+i) because the stress increment is sufficiently small during this interval. 
Nevertheless, the weight, w, would have no significant bearing on the damage 
computation if the difference between dSe/dt and dlWdt derived in eqs. (4- 14a) and (4- 
14b) is small. 
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If a change in the reference point occurs at x r e (tk, tk+i), i.e., if a small cycle is 
closed at x r , then the integral in eq. (4.17) can be split into two parts to exclude the 
point at x r where a step change in c r occurs. 


= f r 
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w — -+(l-w) — - 
dt dt 
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w — - + (1- w) — — 
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dt 


= w ( 8e(a(X r ), Ctrl) - 5e(a(tk), Ori) ) 

+ (1-W) ( 8p(o(X r ), Ori) - 5p(0(tk), Orl) ) 

+ w ( 6 e (a(tk+l), Gr2) - 5e(a(Xr), Gil) ) 

+ (1-w) ( 6p(a(tk+i), Gil) “ 8p(a(x r ), ) 


(4.18) 


where a t i is the reference stress associated with the small cycle, and g t i is the 
reference stress associated with the large cycle.. Figure 4.4 shows both stress-time and 
stress-strain curves when a small cycle is closed between tk and tk+i. If more than one 
cycle is closed during the interval (tk, tk+i), the integration can be split into as many 
parts as necessary to exclude those points where the cycles are closed because these 
points form a set of zero Lebesgue measure as discussed earlier. 



Figure 4.4 Integration When a Cycle Is Closed Between tk and tk+i 
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4.3. Results of Simulation Experiments 

The continuous-time damage model was verified by using the experimental data 
generated in the SAE cumulative fatigue damage test program [Tucker and Bussa 
(1977)]. This program was conducted using a notched specimen shown in Figure 4.5, 
and the materials under tests were Man-Ten and RQC-100 steels which are commonly 
used in the ground vehicle industry. Table 4.1 lists the fatigue properties for these two 
materials. The profiles of load spectrum were recorded from the components operating 
under the actual service condition. Three load profiles, named as suspension, bracket, 
and transmission, were obtained with compressive, approximately zero, and tensile 
mean stresses, respectively. Each load profile was normalized such that the absolute 
value of the maximum load is equal to 999, and small cycles were filtered out of the 
original history. The notch root strain versus applied load was measured from the strain 
gauge, and the relationship is plotted in Figure 4.6. The applied load is converted into 
local strain following Figure 4.6 and the local stress is calculated from the cyclic stress- 
strain curve in eq. (4.4). The continuous-time damage model developed in section 4.2 
is then used to predict the accumulated damage and damage rate. Three different load 
levels of die transmission history are simulated based on the Man-Ten steel data, and 
the results are shown in Figure 4.7, 4.8 and 4.9. The damages computed from the 
strain-life approach are also included for comparison. In the results obtained by the 
strain-life approach, only those points on the peak of a cycle are evaluated, and the 
damages of reversals are computed whenever cycles are not closed. 

Table 4. 1 Smooth Specimen Stress-Strain and Fatigue Properties 


Property 

Man - Ten 

RQC - 100 

Modulus of Elasticity, E 
Cyclic Strain Hardening Exponent, n' 
Cyclic Strength Coefficient, K' 
Fatigue Strength Coefficient, o’ 
Fatigue Strength Exponent, b 
Fatigue Ductility Coefficient, e' 
Fatigue Ductility Exponent, c 

203,000 Mpa 
0.193 
1190 Mpa 
930 Mpa 
-0.095 
0.26 
-0.47 

203,000 Mpa 
0.100 
1150 Mpa 
1165 Mpa 
-0.075 
1.06 
-0.75 


As seen in Figure 4.7, 4.8 and 4.9, the predicted damage generated from the 
continuous-time damage model agrees closely with that from the strain-life approach. 
Each figure shows the accumulated damage of one block load history which contains 
1708 reversals. For high cycle fatigue, Le., maximum load = 15.6 KN, the continuous- 
time damage model tends to predict slightly higher damage than the strain-life 
approach. For low cycle fatigue, i.e., maximum load - 71.2 KN, the continuous-time 
damage model predicts less damage. The relative error, which is defined as the 
difference of damage increments between the two approaches divided by the damage 
increment from the strain-life approach, is within 10% for almost all the reversals in the 
load histories for the two cases where the maximum loads are equal to 71.2 KN and 
35.6 KN. For the maximum load of 15.6 KN, the error is higher since the damage 
computation is more sensitive to a small change of stress or strain in the case of high 
cycle fatigue but the mean value is nearly zero. The results from both approaches are 
considered to be in fair agreement in view of scattering of the fatigue test data. 












24 



Time (Normalized) 



Time (Normalized) 

Figure 4.7 A Comparison of Continuous-Time Damage Model with 
the Strain Life Approach for Maximum Load = 15.6 KN 
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Figure 4.8 A Comparison of Continuous-Time Damage Model with 
the Strain Life Approach for Maximum Load = 35.6 KN 
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Figure 4.9 A Comparison of Continuous-Time Damage Model with 
the Strain Life Approach fen* Maximum Load = 71.2 KN 
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4.4. Modeling of Nonlinear Cumulative Damage Using Damage Curve Approach 

The continuous-time damage model developed above is based on the linear 
damage accumulation following the Palmgren-Miner's rule. Although this concept of 
linear damage accumulation has been widely used due to its simplicity in computation, 
the cumulative damage behavior is actually nonlinear [Suresh (1991)]. Experimental 
results show that the linear damage rule predicts less damage if a few cycles of high 
stress are applied before testing with low stress. This phenomenon is known as the 
sequence effect. In many engineering applications, the components are usually 
subjected to loading with varying amplitudes. Due to this sequence effect, the linear 
rule of damage accumulation, which is commonly used for fatigue life assessment, 
could lead to erroneous results. A nonlinear cumulative damage representation needs to 
be established for more accurate life prediction of the critical components. 

Hie concept of a nonlinear damage curve to represent the damage was first 
conceived by Marco and Starkey (1954). No mathematical representation of a damage 
curve was proposed at that time because the physical process of damage accumulation 
was not adequately understood. Manson and Halford (1981) proposed die double linear 
rule primarily based on the damage curve approach for treating cumulative fatigue 
damage. In their paper, an effort was made to mathematically represent the damage 
curve and approximate it by two piecewise line segments. The total fatigue life was 
thereby divided into two phases so that the linear damage rule could be applied in each 
phase of the life. A concept similar to the damage curve approach was proposed by 
Bolotin (1989) with a mathematical representation which does not necessarily assess 
the damage on the basis of cycles and is more appropriate for modeling in the 
continuous-time setting. Therefore, Bolotin's approach is adopted in this research for 
the development of a continuous-time model with nonlinear damage accumulation. A 
review of the damage curve approach following Bolotin's approach is presented below. 

Figure 4. 10 shows a comparison of the accumulations of linear damage, D/, and 
nonlinear damage, D, as a function of the cycle ratio, n/N, where n is the actual number 
of cycles undergone and N is the number of cycles to failure under a constant amplitude 
load. The damage accumulates along the curve as the loading cycles are applied. For 
example, if a specimen is subjected to ni cycles of a constant stress amplitude, for 
which the fatigue life is Ni cycles, the accumulated damage will be D a indicated in 
Figure 4.10, where the abscissa is the normalized cycle ratio with respect to its fatigue 
life and the ordinate is the damage accumulation of the specimen. Bolotin used the 
following analytical relationship between D and D/: 

D = (D/) Y( ° a) (4.19) 

where the exponent 7 describes the nonlinearity of the curve and usually is a function of 
the stress amplitude c a . 

Eq. (4.19) indicates that the nonlinear damage may accumulate along different 
curves under different stress amplitudes. Generally speaking, for high-strength 
materials that usually strain soften [Hertzberg (1989)], a smaller load amplitude tends to 
make the damage curve more nonlinear, i.e., increase the Y-parameter. To realize the 
effects of nonlinear damage accumulation, consider a smooth specimen be subjected to 
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two levels of cyclic loading, Aai and Aa2 (say Aaj>Ao2), with respective fatigue lives 
of Ni and N2 cycles (Ni<N 2). If Aoi is applied for ni cycles followed by A02 for n2 
cycles when the specimen fails, then the linear rule would underestimate the damage, 
i.e., (nj/Ni + n2/N2)<l. Next consider the reverse situation where A02 is applied first 
for n2' cycles followed by Aai for ni' cycles to failure of the specimen. The linear rule, 
in this case, overestimates the damage, i.e., (ni’/Ni + n2'/N2>>l. Suppose OXZ and 
OYZ are the nonlinear damage curves for cyclic loads, Aai and Aa2 respectively, as 
shown in Figure 4 . 11 . First, let ni cycles of Aai be applied prior to A02. As the load 
cycles are applied, the accumulated damage follows the curve OXZ and reaches the 
point A with its value equal to D a . Then, let A02 be applied until the component fails. 
The damage accumulation continues from its present state D a along the curve OYZ 
until the nonlinear damage D reaches 1 . That is, the damage accumulates from point B 
to Z. Suppose the number of cycles needed for the component to fail under A02 is n2. 
As indicated in Figure 4 . 11 , the linear rule underestimates the damage by the length 
AB, which is equal to l-ni/Ni-n2/N2. In the reverse case, the linear rule overestimates 
the damage by the length of AH' as shown in Figure 4 . 12 . 



Cycle ratio n/N 


Figure 4.10 Nonlinear Damage Curve under Constant Amplitude Loading 


The above example illustrates the basic concept of the damage curve approach 
for two-level loading. In the case of multiple-level loading, the same procedure can be 
applied by identifying the present damage state and following the damage curve 
associated with the current loading condition. If a component is subjected to spectral 
loading, then the techniques of cycle-counting, prediction of linear damage increments 
and computation of nonlinear damage via the damage curve approach need to be 
integrated into a single procedure. 
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In the damage curve approach described above, the Y-parameter is usually assumed to 
be dependent solely on the stress amplitude level [Manson and Halford (1981)]. High- 
strength materials such as 4340 steel usually yield very large values of Y especially 
under high-cycle fatigue. Manson (1981) made an attempt to fit the Y-parameter in a 
damage model that is structurally similar to eq. (4.19) and includes an initial damage 
term. 



Linear Damage Accumulation D/ 
Figure 4.11 Nonlinear Damage Accumulation (Aai > A 02 ) 



Linear Damage Accumulation D; 


Figure 4.12 Nonlinear Damage Accumulation (A 01 < A 02 ) 
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The y-parameter was found to have the following form: Y = (2/3)Nf°- 4 , where Nf 
is the number of cycles to failure under a constant load amplitude. For example, if Nf is 
assumed to be 10 5 cycles, the Y-parameter is then equal to 66.7 which appears to be 
extremely high especially for damage computation when D«l. It also follows from 
eq. (4.19) that a large Y, at an early stage of fatigue life, shall yield a very small damage 
which could be out of the range of the computer precision. This necessitates die 
formulation of a computationally practical method of damage prediction. 

The cause of large Y-parameter is that, in the conventional damage curve 
representation as shown in eq. (4.19), the values of Y are assumed to be constant at 
certain stress amplitudes. These curves do not accurately describe the nonlinearity of 
damage accumulation process for the entire fatigue life of the component It follows 
from a crack propagation model such as the Paris model [Paris and Erdogan (1963)] 
that the crack growth rate is dependent not only on the stress amplitude but also on the 
current crack length. Recognizing the fact that the crack itself is an index of 
accumulated damage, it is reasonable to assume the Y-parameter to be dependent on 
both stress amplitude and the current level of damage accumulation., i.e., Y - Y(a a , D). 
Therefore, eq. (4.19) should be modified as 

D = (D / ) Y(a *’ D) ( 420 ) 

where D and D/ are the current states of nonlinear and linear damage accumulation, 
respectively. Although the above eq. (4.20) has an implicit structure, it can be solved 
via a recursive relationship. 

The next part of this section describes a modification of the damage curve 
approach to develop a nonlinear damage model in the continuous-time setting. 
Following the concept of the linear damage model in the continuous-time setting, the 
nonlinear damage at any point on a rising reversal can be obtained as explained below. 

Referring to the bottom part of Figure 4.13, let A be any point on the rising 
reversal and R be its reference point as determined from the rainflow cycle counting 
method [Dowling (1983)]. Let the current state of damage at the reference point R be 
equal to D r and let ORAZ be the damage curve associated with the stress amplitude 
RA/2 as shown in the top part of Figure 4.13. Corresponding to the nonlinear damage 
D r , the notion of the "virtual" linear damage, D&, is brought in as follows: 

l 

D, a = (D r ) 7 ‘ («1) 

where Yr * s the Y-parameter associated with the stress amplitude RA/2 and nonlinear 
accumulated damage Dr at point R. The tom "virtual" means that D/ a is the linear 
damage which would be incurred if the component had been subjected to the cyclic 
stress of constant amplitude, RA/2, from its initial damage state to the current damage 

state at R. Similar to Yr “ eq* (4.21), the Y-parameter associated with the damage state 
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D a at the point A is defined as yf. Referring to Figure 4.13, the functional 

relationships among D a , Dr and the corresponding Y-parameters y a and y? are defined 
as follows: 


Yl - Y(®a» D a) [y? = Y(o a ,D r ) 

l D a=(D/a+AD/ a ) r * Pr = (D/a) ' 


(4.22) 


where o a is the stress amplitude RA/2, and AD/ a is the linear damage increment 
between R and A obtained via the procedure described in Section 4.2. The damage Dr 

Z 




Strain 


Figure 4.13 Nonlinear Damage Increments 
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in the right part of eq. (4.22) represents the damage state at the reference point R and 
therefore its value is already computed from the past load history. Having known Dr, 

the Y-parameter yf and the "virtual" linear damage D/ a at R for the stress amplitude 
RA/2 can be evaluated from the right hand part of eq. (4.22). The two unknowns, D a 

and y|, which represent the damage state and the Y-parameter at the current point A, are 

computed by solving the equation pair in the left part of eq. (4.22). Now, the 
(nonlinear) damage increment from the point R to A in Figure 4.13 can be computed as: 

AD a =D a -D r = (D <a + AD <a ) r * - D <a tf (4.23) 

In summary, the accumulated damage at any point within a reversal can be 
obtained by solving the following nonlinear equations: 

Y = Y(c?a> D ) (4.24) 

D = (D/ + AD/)* (4.25) 

where D/ is the "virtual" linear damage at the reference point and AD/ is the linear 
damage increment for the stress amplitude o a relative to the reference point 

The pair of equations (4.24) and (4.25) does not have a closed form solution and 
therefore needs to be solved by an iterative method. However, an iterative solution at 
every point in the load history may not be practical from the perspective of 
computational efficiency. An efficient approach of obtaining a numerical solution 
would operate on a set of discrete points in the load history such that there is only a 
very small increment of damage between two consecutive points. Thus, both eqs. 
(4.24) and (4.25) can be treated as linear line segments between the points R and A. 
This assumption is valid during the entire loading history with the possible exception of 
very low-cycle fatigue. One more interesting observation is that the Y-parameter is 
generally a monotonically increasing function of the nonlinear damage D. This can be 
interpreted from the Paris equation that the growth rate of a macrocrack becomes larger 
as the crack length increases. Therefore, the damage rate would be larger at a higher 
degree of nonlinearity, which implies a larger value of Y. This phenomenon, however, 
may not be true if the stress intensity factor range is below the long crack threshold or if 
the material strain hardens. As seen in Figure 4.14, for the general case of Y>1 (e.g., 
high strength materials that usually strain soften [Hertzberg (1989)]), Y is a 
monotonically increasing function of D in eq. (4.24), and D is a monotonically 
decreasing function of Y in eq. (4.25). On the other hand, if Y is less than 1 (e.g., ductile 
materials that usually strain harden), the characteristics of both eqs. (4.24) and (4.25) 
could be reversed. 

Having computed the linear damage D/ and linear damage increment AD/, a 
procedure for solving eqs. (4.24) and (4.25) to obtain the nonlinear damage, D, is 
described below: 


1. LetDi = D r and Yi=Y r . 

2. Find D2 = (D/4AD/ )Yi from eq. (4.25) and Y2 = Y (o a , P2) from eq. (4.24) 
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3. Compute Y3 as Y3 = Yi x — - — - — - from eqs. (4.20) and (4.24) 

l°g(D/ + AD^) 

4. Find D a in Figure 4.14, which is approximated as the point of intersection of two 
straight line segments: 


D Di(Y 2 -Yi)+D 2 (Y 3 -Yi) 
a Y 3 +Y 2 -2Y! 


(4.26) 


The above procedure computes the nonlinear accumulated damage at any point 
on the rising reversal. If the stress is monotonically decreasing, there is no damage 



Figure 4.14 Computation of Nonlinear Cumulative Damage 

increment as described in eqs. (4.14) and (4.15) for the case of linear damage 
accumulation. Finally, the rate of nonlinear damage accumulation is obtained directly 
by differentiating eq. (4.20) with respect to time t 


?-*>« r 



(D^lnD,*^. 


(4.27) 


4.5. The Y -Parameter Fitting for the Nonlinear Cumulative Damage Model 

One major task in the above approach is to identify a mathematical 
representation for the Y-parameter as a function of the applied stress amplitude and the 
current damage state. It requires the knowledge of physical process of damage 
accumulation which may be obtained from either experimental data or a combination of 
experimentation and analysis with an appropriate definition of damage. The Y- 
parameters have different values, in general, for different materials and follow different 
equation structures. Furthermore, because the mechanisms attributed to the damage 
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accumulation at various stages of fatigue life are different, no single approach can 
provide accurate assessment of damage throughout the fatigue life of a component. It is 
difficult, if not impossible, to construct a single structure for representation of Y. An 
alternative approach is to evaluate Yby interpolation. Once the damage is appropriately 
defined and an analytical method is selected, the information needed for the nonlinear 
damage model can be generated via experimentation or analysis for constant stress 
amplitudes. The values of Y are then computed by eq. (4.20) for a given stress level and 
the damage data ranging from Do (initial damage state) to die failure condition at D=l. 
These generated data, Y versus D, are then plotted for various amplitudes of stresses. 
These curves can be either fitted by nonlinear equations, if possible, or as linear 
piecewise representations. Since it is not practical to perform experiments at 
infinitesimally small increments of stress amplitude, the values of Y can be only 
experimentally determined at selected discrete levels of stress amplitude. The values of 
Y for other stress amplitudes can then be interpolated. Since the characteristics of Y may 
strongly depend on the type of the material, availability of pertinent experimental data 
for the correct material is essential for the damage-mitigating control synthesis. The 
remaining part of this section presents the results of Y-parameter fitting based on the 
experimental data [Swain et al. (1990)] for the material AISI 4340 steel. 

In the model of Newman et al. (1981), the stress intensity range, AK, used in the 
Paris equation is replaced by effective stress intensity range AKeff. The crack opening 
stress is determined by a crack closure model which is similar to the Dugdale model 
[Dugdale (I960)] but it is modified to leave plastically deformed material in the wake 
of the advancing crack tip. In this simulation, however, the crack opening stress is 
computed by simplified equations [Newman (1984)] which are obtained through curve 
fitting based on the original model. The experiments show that a unified approach 
[Newman et al. (1992)] based on the crack closure concept can be used for damage 
prediction starting from the initial defect (microcrack) to the failure of materials 
without significant errors. A relationship between the effective stress intensity factor 
range and crack growth rate obtained from the experimental data of AISI 4340 is given 
in Table 4.2 [Swain et al. (1990)]. 

Table 42 Effective Stress Intensity Factor Range 
versus Crack Growth Rate Relationship 


AKeff 

MpaVm 

da/dN 

m/cycle 

3.75 

3.0E-10 

5.30 

2.0E-09 

7.30 

7.0E-09 

15.00 

4.5E-08 

50.00 

5.5E-07 

120.00 

3.0E-05 


In Table 4.2, the values of AKeff and da/dN are linearly interpolated between 
two consecutive data points in the logarithmic scale while those beyond the maximum 
and minimum are extrapolated. In the simulation, the fatigue life is first predicted using 
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the strain-life approach. This fatigue life is assumed to be identical to that obtained via 
the crack closure model, and is used to identify the initial crack size. The damage, in 
this case, is defined as the normalized crack length with respect to its critical length to 
failure (i.e., D=a/a*) where the critical crack length a* is obtained from fracture 
toughness of the material. The material and fatigue properties of AISI 4340 steel are 
given in Table 4.3 [Boiler and Seeger (1987)]. 


Table 4.3 Material Properties of AISI 4340 


Young's Modulus (E) 

193500 N/mm^ 

Yield Strength, Monotonic (Gy) 

1374 N/mm2 

Yield Strength, Cyclic (Gy' ) 

905 N/mm^ 

Cyclic-Strength Coefficient (K*) 

1890N/mm2 

Cyclic-Strength Hardening Exponent (n’) 

0.118 

Fatigue-Strength Coefficient (of 1 ) 

1880 N/mm2. 

Fatigue-Strength Exponent (b) 

-0.086 

Fatigue-Ductility Coefficient (Ef*) 

0.706 

Fatigue-Ductility Exponent (c) 

-0.662 


The results of the nonlinear damage D versus the linear damage D/ which is 
essentially the cycle ratio are used to compute Y from eq. (4.20). Figure 4.15 shows the 
relationship between Y and D for different values of the (constant) stress amplitude G a 
as a series of curves in the logarithmic scale. If these curves are generated to be closely 
spaced, the values between two curves can be obtained via linear interpolation without 
any significant error. Once the Y-parameters are computed as a function of two 
independent variables o a and D, the nonlinear damage model described above can be 
readily used to simulate the nonlinear damage behavior. The nonlinear damage 
accumulations corresponding to different stress amplitudes are plotted in Figure 4.16. 

It is seen in Figure 4.15 that Y-parameter is strongly dependent on the current 
damage level in contrast to the conventional damage curve approach where Y i s 
assumed to be constant At the initial stage of fatigue life, Y-parameters are only 
modestly larger than one so that it is more realistic for implementation in computer 
simulation. In general, the values of Y are larger for smaller stress amplitudes. This 
indicates higher degrees of nonlinearity in high cycle fatigue as can be seen from the 
damage curves in Figure 4.16. The accumulated damage is much smaller in the early 
stage of fatigue life compared to the linear damage accumulation and accelerates much 
faster when the damage is close to failure. This implies a significant improvement in 
fatigue life prediction based on the nonlinear damage rule if the Y-parameters are 
properly selected from the experimental observations or on the basis of analytical 
results. 
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4.6. Continuous-Time Damage Modeling of Fatigue Crack Growth 

The concept of a continuous-time damage model, developed in Section 4.2 
using the strain-life approach, can be extended to formulate a continuous-time model 
for crack propagation. In the most commonly used fatigue crack growth models, the 
growth rate of the crack length, a, relative to the number of cycles, N, is directly 
dependent on the stress intensity factor range, AK. Most of these models can be 
expressed in a mathematical form as: 

^ = f,(AK) = f(Ao,a) (4.28) 

dN 

where AK is a function of both stress amplitude and crack length. In this report, the 
crack closure model is adopted and therefore AK is replaced by the effective stress 

intensity factor range AKeff and Ao=o-o 0 , where o 0 is the crack opening stress. To be 
consistent with the strain-life approach, damage is assumed to be identically equal to 1 
when the crack length reaches the critical value, a*, which can be estimated from the 
fracture toughness of the component The damage is assumed to occur only if the 
instantaneous stress is increasing and the crack opening stress, o 0 , is exceeded 
[Newman, 1981]. This implies that the damage rate is equal to zero during unloading. 
The damage increment AD during the rising reversal can be expressed as a function of 
Aa and D: 



-L x — = -i- f (Aa,a*D) if octo 
a dN a 


(4.29) 


where the fatigue damage is defined as D=a/a* in the absence of a more precise 
definition of damage. Since the damage increment is relatively small, the total damage 
D is assumed to remain constant during a reversal. The damage rate can be obtained by 
differentiating eq. (4.29) with respect to time. 


— = 4-x df(A °’ a * D) x— (4.30) 

dt a do dt 

A simple example based on the Paris model [Paris and Erdogan (1963)] is given below. 

Let [tk] be the sequence of time instants when damage values are to be 
estimated, and (a(tk)} be the corresponding sequence of stresses. Then, the Paris 
model can be modified to yield the time derivative of dama ge as: 


dD C / \n-lf I *77—/ * \Y*do 
— = -^(a-o 0 ) [Vita D F(a D,w)J — 


a 

®=0 

dt 


if o(tk) xr(tk-i) and o(tk) > <J 0 


otherwise 


(4.31) 
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where F is the correction factor depending on the geometry and w is the width of the 
specimen. The application of eq. (4.31) requires estimation of the crack opening stress, 
a 0 , which can be calculated from the crack closure model [Newman (1981)] based on 
the load history. Following a simplified method [Newman (1984)] for computing ao, 
the damage increment between two consecutive instants tk-i and 4 can be obtained by 
integrating eq. (4.31) as follows: 

= 4[(o<t k ) - o 0 )" - (Max(c(t k _i ), o 0 ) - o 0 )“] x (Vjca*D(t k _i ) F(a*D(t k _iXw))" 

if o(tk) xjfe.i) and o(tk) > a 0 

D(t*) = D(tk-i) otherwise. (4.32) 



CHAPTER 5 


SIMULATION RESULTS AND DISCUSSIONS 

The damage mitigation concept, described in Chapter 2, has been verified by 
simulation experiments for open loop control of a reusable rocket propulsion engine 
such as one described in [Sutton (1992); Duyar et al. (1991)]. The plant model under 
control is a simplified representation of the dynamic characteristics of a bipropellant 
rocket engine as shown schematically in Figure 5.1. The prebumer serves as a gas 
generator for driving the liquid hydrogen (LH 2 )-fuel turbopump. In this model, oxidant 
is separately supplied to the prebumer and the main combustor chamber. Standard 
lumped parameter methods have been used to model the nonlinear plant dynamics in 
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Figure 5.1 Schematic Diagram of a Bipropellant Rocket Engine 
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the state-space form where the plant state vector consists of turbine shaft speed, pump 
(LH 2 -fuel) mass flow rate, prebum er gas pressure, prebumer gas density, combustor 
gas pressure, combustor gas density, and the two flow rates of oxidant into the 
prebumer and the main combustor, respectively. The critical plant output variables are 
combustor gas pressure and the oxidant/ fuel (O 2 /H 2 ) ratio, and the control inputs are 
the areas of the two oxidant valves. The governing equations for the lumped parameter 
model of the plant are given in Appendix C. The structural model, as delineated in 
Appendix D, calculates the cyclic stresses at the root of a typical turbine blade that is 
presumed to be a critical component The blade is represented by a three-node beam 
model with six degrees of freedom at each node while the first node is kept fixed. The 
load on the blade is assumed to consist of two components. The first component is due 
to the (time-dependent) drive torque which is derived as an output of the plant model. 
The second component is a dynamic term which represents the oscillatory load on the 
blade as it passes each stator. It is the second component that causes high cycle fatigue 
at the root of the blade while the first component is largely responsible for the mean 
stress. The fatigue damage model formulated in Chapter 4 was used to generate the 
results in Figures 5.2 to 5.14. 

5.1. Simulation of the Bipropellant Rocket Engine under Open-Loop Control 

The simulation experiments serve to evaluate the plant dynamic performance 
and damage of the critical component when the oxidant valves are manipulated to vary 
the engine thrust following the open-loop control policy developed in Chapter 3. To 
demonstrate the broad concepts of fatigue damage mitigation, the nominal plant model 
was used in the simulation experiments with exact initial conditions and no 
disturbances and noise. However, if these conditions are not met, additional feedback 
control will be necessary because the open-loop control alone would be inadequate for 
plant operations as discussed in Appendix A. Following the structure in eq. (3.3), the 
cost functional J for nonlinear programming was selected to generate the open loop 
control policy as: 


J = 


N-l 

£ [x k T Qx k +v k T Sv k +u k T Ru k + W(0 2 /H 2 ) k 2 ] 
k=0 


(5.1) 


where the deviations, x k and u k , in the plant state vector and the control input vector 
are as defined in eq. (3.3); and the diagonal matrices Q, S, R and the scalar W serve as 
relative weights of the individual variables. Since the rocket engine performance is 
very sensitive to the oxidant/fuel (O 2 /H 2 ) ratio, it was brought into the cost functional in 
eq. (5.1) to prevent any large deviations from the desired value through the transients. 
If the main combustion chamber is selected as one of the critical components for 
damage mitigation, then the gas temperature has to be controlled within a small bound. 
In that case, the cost functional in eq. (5.1) need not be explicitly dependent on the 
O 2 /H 2 ratio which is directly related to the main chamber temperature in the range of 
normal operations. Simulation results were obtained under the following conditions: 

Chamber pressure weight Qs5=12; and all other weights Qii=l, i*5; 

Control input weight R=I where I is the identity matrix; 

Damage Rate Weight Sa=0, i=l, *•*, N; and O 2 /H 2 ratio weight W=10. 
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The rocket engine model is initiated from an initial equilibrium condition at 
2700 psi chamber pressure and O 2 /H 2 ratio of 6.02. From this condition the 
optimization procedure steers the plant to a new equilibrium position at 3000 psi and 
the same O 2 /H 2 ratio of 6.02 in 50 milliseconds. The control commands to the two 
oxidant valves, are updated at every one millisecond. That is, N is equal to 50 in eq. 
(5.1). The performance cost to be minimiz ed is based on the deviations from the final 
equilibrium condition at 3000 psi. The results of simulation experiments for the two 
following conditions are presented as series of curves in Figures 5.2 to 5.14: 

Simulation Condition 1 : Initial damage is set to Do=0.01 for the unconstrained case and 
two constrained cases in which the damage rate constraints are listed in Table 5.1. 

Simulation Condition 2 : To examine the effects of initial damage on nonlinear damage 
accumulation, simulation results are generated under three different initial values of the 
accumulated damage, namely, Do=0.001, 0.01 and 0.1 while the damage rate constraint 
is set identical to that of the constrained case 1 in Table 5.1. 

For both conditions 1 and 2, the accumulated damage (<p) was not constrained. 

Table 5.1 The Constraints for the Two Cases under Simulation Condition 1 


Time 

0 ms to 2 ms 

2 ms to 3 ms 

3 ms to 50 ms 


Case 1 

Damage rate constraint (P) 
l.OxlO* 6 sec' 1 
2.5x10-6 sec 1 
5.0x10*6 sec- 1 


Case 2 

Damage rate constraint (P) 
0.2x10-6 sec 1 
0.5x10-6 sec 1 
l.OxlO- 6 sec 1 


52 . Results and Discussions 

The transient responses in Figures 5.2 to 5.14 examine the various engine 
variables under the above two simulation conditions. Figures 5.2 and 5.3 show the 
transient responses of the two oxygen valves resulting from the optimization over the 
time frame of 0 to 50 ms where the control action is updated at every millisecond. The 
corresponding changes in the oxidant flows are seen in Figures 5.4 and 5.5. Figure 5.6 
also shows the transients of the fuel flow (Le., liquid hydrogen). The responses of both 
control valves (and therefore the oxidant flows) become more restricted and sluggish as 
the damage rate constraint is made stronger. Similar effects are observed when the 
initial damage is increased. The rationale for this behavior is that, for a given stress 
amplitude, the damage rate increases with an increase in the initial damage. This 
dependence on the initial damage results from the Y-parameter in the nonlinear damage 
model as defined in eq. (4.19) in Chapter 4, and does not occur in the linear damage 
model where Y is identically equal to 1. It is important to note that Y is greater than 1 in 
this study because the AISI 4340 steel used here is a high strength material. However, 
for ductile materials such as copper-base alloys used in the thrust chamber cooling 
tubes, Y may be less than 1 due to strain hardening, and therefore the dependence of the 
control commands on the initial damage could be significantly different 
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Figures 5.7 to 5.11 exhibit the effects of the varying oxygen inlet flow to the 
prebumer and the main combustor on the engine dynamics. The resulting transients of 
the process variables, namely, O 2 /H 2 ratio, and die pressure and temperature in the 
prebumer and main combustion chamber, are shown for the two simulation conditions. 
As expected, for a given level of initial damage, both pressure and temperature 
dynamics tend to be slower as the constraint is made more severe. Similar effects are 
seen by increasing the initial value. Do, of the accumulated damage. The combustor 
pressure is seen to rise monotonically in all cases after a small dip at about 2.5 ms while 
the prebumer pressure keeps on increasing. These plots are largely similar except for 
the transients from 1 ms to 10 ms. Virtually all of fatigue damage accumulation in this 
transient operation takes place during this short interval as seen in Figure 5.14. 
Furthermore, the net excursion of the O 2 /H 2 ratio is in the range of 6.0 to 6.5 in all 
cases for the up-thrust transient of the rocket engine. The 6.5 mixture ratio is about the 
limit that would be tolerated during a transient excursion. The overshoot in the mixture 
ratio is caused by a drop in the hydrogen flow as seen in Figures 5.6 and 5.7. At this 
point the turbopump demands more torque to increase its speed so that the pump 
pressure can be elevated to generate a higher value of hydrogen flow for the desired 
mixture ratio. This results in a peak overshoot in the mean stress as shown in Figure 
5.12. This sharp increase in stress is the major cause of enhanced damage in the turbine 
blades. 


Since the turbine blades are the critical components for damage analysis in this 
study, the pressure ratio across the turbine which directly influences the torque is very 
important. For a given prebumer pressure, as shown in Figure 5,9, a reduction in the 
combustor pressure causes an increase in the turbine torque. It is the turbine torque and 
speed that set the stress and fatigue damage on the turbine blades. Therefore, the dip in 
the combustor pressure at about 2.5 ms is largely responsible for the peak mean stress 
displayed in Figure 5.12. Both the mean stress and stress amplitude are the basic inputs 
to tiie damage model. 

The graphs in Figures 5.12 to 5.14 compare the damage rate and the 
accumulated damage for the two simulation conditions along with the transient 
responses of the mean stress. For the unconstrained case, the peak stress causes the 
largest overshoot in the damage rate which is plotted on a logarithmic scale. In 
contrast, for the initial damage of 0.001, the damage rate is within the limit of the 
constraint even though the peak of mean stress is the largest This phenomenon is a 
result of a relatively small slope of the damage curve at the initial stages of the fatigue 
life. The accumulated damage, plotted on a linear scale, is seen to be significantly 
influenced by the constraints and also by the initial damage. This suggests that for 
reusable rocket engines, the constraints need to be appropriately specified based on the 
knowledge of the initial damage. The damage rate is dependent on the sequences of 
control commands, and the oxidant flows into the prebumer and combustor are changed 
in response. Therefore, if the initial damage cannot be accurately assessed, then it 
might be safe to generate the control command sequences on the assumption of a 
conservative, i.e., larger, value of the initial damage at the expense of the engine 
performance. 

The important observation in these simulation experiments is the substantial 
reduction in the accumulated damage, thereby extending the service life of the turbo- 
pump. The accumulated damage in the unconstrained case is seen to be about four to 
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twelve times that of the constrained case. This is a clear message that the consideration 
of damage in the control of transients to which a rocket engine is exposed can have a 
considerable impact on the life of critical components (in this case, the turbine blades). 
It is noted that there is practically no penalty in the response times of chamber pressure, 
i.e., the engine thrust, between the unconstrained and the constrained cases. If one is 
willing to pay a small price in response time, much larger gains on damage 
accumulation can be achieved. 

Figures 5.2 to 5.14 exemplify the effects of upthrust transients during a short 
period of 50 ms. Complete operations of a rocket engine during a single flight include 
many such upthrust transients, and the steady-state operation may last for several 
hundreds of seconds. Although the damage rate during the steady state is much smaller 
than that during a transient operation, the total damage accumulation during the steady 
state may not be relatively insignificant Therefore, during a flight of the (reusable) 
rocket engine, the cumulative effects of the transient and steady state operations need to 
be considered in the optimization procedure as discussed in Chapter 3. 

The simulation experiments, described above, only consider a single point of 
critical stress, namely, the turbine blades. In this case, the damage vector is one- 
dimensional. Simultaneous control of damage at several other critical areas in the 
rocket engine, such as the nozzle lining, shall render the damage vector to be multi- 
dimensional. The optimization problem is then to generate a control sequence that will 
not only make a trade-off between the performance and damage but also strike a 
balance between potentially conflicting requirements of damage mitigation at the 
individual critical points. 
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Figure 5.3 Transient Responses of Combustor O 2 Valve Areas 
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Figure 5.4 Transient Responses of Prebumer O5 Inlet Flow 
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Figure 5.5 Transient Responses of Combustor 0 2 Inlet Flow 
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Figure 5.6 Transient Responses of H 2 Mass Flow 
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Figure 5.7 Transient Responses of O 2 /H 2 Ratio 
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Figure 5.1 1 Transient Responses of Combustor Temperature 
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CHAPTER 6 


SUMMARY, CONCLUSIONS AND RECOMMENDATIONS 
FOR FUTURE RESEARCH 


The theme of the research in damage-mitigating control of mechanical systems 
is summarized as follows: 

• For control of mechanical systems, damage prediction and damage mitigation arc 
carried out based on the available sensory and operational information such that the 
plant can be inexpensively maintained, and safely and efficiently steered under 
diverse operating conditions. 

• High performance is retained without overstraining the mechanical structures such 
that the functional life of critical components is increased resulting in enhanced 
safety, operational reliability, and availability. 

In this report, a fatigue damage model has been developed in the continuous- 
time setting such that this damage model can be integrated with the plant dynamic 
model for die purpose of control synthesis. A procedure of synthesizing an optimal 
policy for open-loop control has been formulated via nonlinear programming under the 
constraints of damage rate and accumulated damage. Finally, efficacy of the proposed 
damage-mitigating control system has been demonstrated via simulation of the transient 
operations of a reusable rocket engine. 

This research in damage-mitigating control is interdisciplinary and addresses the 
fields of active control technology and structural integrity. Extended service life 
coupled with enhanced safety and high performance will have a significant economic 
impact in diverse industrial applications. Examples include reusable rocket engines for 
space propulsion, rotating and fixed wing aircraft, fossil and nuclear plants for electric 
power generation, automotive and truck engine/transmission systems, and large rolling 
mills. Furthermore, as the science and technology of materials evolve, the updated 
damage characteristics of the structural components can be systematically incorporated 
within the framework of the proposed control system. 


6.1. Recommendations for Future Research 

The structure of the proposed damage-mitigating control system is built upon 
the concepts of both the conventional state feedback control systems and fatigue life 
prediction. Certain assumptions and approximations are made in the modeling of 
damage dynamics for practicality of implementation. Further research is necessary for 
more accurate damage prediction of the critical components. In the synthesis of robust 
damage-mitigating control, research must be conducted for accommodation of the 
disturbances and errors in modeling of both plant dynamics and damage. The following 
topics of research are recommended in the areas of damage prediction and control. 
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1. Identification of material parameters an d verification of the damage model 

The continuous-time damage model is built upon the strain-life approach which 
assumes that the coefficients for material properties are constant and cyclic stress-strain 
behavior remains stabilized for most of the fatigue life. The computation load for real- 
time damage control could be reduced by assuming the parameters of the structural 
materials to be invariant relative to temperature. However, if a service component is 
exposed to high temperature, then the effects of parametric variations and additional 
phenomena of damage such as creep and corrosion must be taken into consideration. 

Apparently no consistent definition of damage exists in the current literature. 
Therefore, damage must be precisely defined for accurate prediction of service life 
especially if the synergistic effects of fatigue and other phenomena such as creep and 
corrosion are considered. In this research, the linear damage for a given stress or strain 
amplitude has been defined on the basis of the total number of cycles to failure. This 
definition has been extended for nonlinear damage where the concepts of normalized 
crack length and crack propagation have been used to identify the Y-parameter of the 
damage curve. Although tins model of damage prediction is considered to be adequate 
for a certain range of service life, it may not be sufficiently accurate for the entire 
range because of the different mechanisms of crack initiation and growth. For an 
appropriate damage definition, experimental research is needed for a better 
understanding of the micro-structure and failure mechanisms of the structural materials. 

2. Creep-fatigue interaction and formulation of a multi-dimensional damage vector 

The damage equations developed in the main text of this report consider only 
the fatigue failure of the materials. If the service components are exposed to high 
temperature or gaseous environment, the creep or corrosion may play an important role 
in the failure of the materials. The damage prediction has to incorporate the effect of 
creep-fatigue interaction for more accurate results. A general structure for modeling of 
creep-fatigue and corrosion-fatigue failure is proposed in Appendix E. If the creep or 
corrosion damage becomes important, or several critical points for failure are 
considered, then the damage vector should be multi-dimensional and needs to be 
optimized to achieve the best trade-off between the damage of individual critical 
components and the system performance. 

3. Stochastic modeling of damage dynamics 

In contrast to the linear damage rule, the nonlinear damage model developed in 
Chapter 4 requires the information of initial damage to initiate the damage computation. 
The initial damage in a specific component, however, is difficult to be identified 
because the size, distribution, and geometry of the defects in the material may vary over 
a wide range. Since a large part of the fatigue life may be consumed during the initial 
stage of so called crack initiation especially for the components with smooth geometry 
or under light load, a small error in the estimation of initial damage may cause a 
significant difference in the fatigue life prediction. Other parametric and non- 
parametric uncertainties in the damage model such as those resulting from variations in 
the material properties also contribute to the inaccuracy and scattering in prediction of 
the service life. If an average value is used in the deterministic fatigue damage model, 
the variance of the predicted life is likely to be large and therefore the damage model 
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would have poor repeatability. However, the modeling uncertainties could be directly 
taken into account if the damage dynamics are modeled by stochastic differential 
equations. This approach could minimiz e the modeling error in a statistical sense and 
be used in fatigue failure prognosis and risk analysis. 

4. flgm b lggLg fflffi fl 

The open-loop control policy optimizes the system performance under the 
damage constraints without considering the plant modeling errors, disturbances, and 
sensor noise. The closed-loop control is required to maintain the system along the 
desired trajectory or at the normal steady-state operating condition. A closed loop 
system could be designed on the basis of the linearized plant model at the normal 
operating point by using the techniques of robust control synthesis such as LQG/LTR, 
H2/H00, or ji-synthesis. However, to ensure that the feedback control does not violate 
the damage constraints, a two-tier structure of the control system might be needed. The 
concept of this hierarchically structured closed-loop control system is outlined in 
Appendix A. 



APPENDIX A 


CLOSED LOOP CONTROL 

Nonlinear progr ammin g generates an open-loop control policy to achieve 
optimal performance under the specified constraints of damage rate and accumulated 
damage as described in Chapter 3. However, because of plant modeling uncertainties 
(including unmodeled dynamics), sensor noise and disturbances, the actual plant 
response shall deviate from that of the modeled system when the plant is excited by the 
sequence of open-loop control commands. Therefore, a closed-loop control system is 
necessary to compensate for these deviations, and a state feedback controller may serve 
this purpose. If all plant states are not measurable or if the sensor data are noise- 
contaminated, a state estimator (e.g., a minimum variance filter) is necessary to obtain 
an estimate of the plant state vector. If the deviations from the nominal trajectory are 
not large, the state feedback and state estimator gain matrices could be synthesized 
based on a linearized model of the plant. However, if the plant is required to be 
operated over a wide range (for example, scheduled shutdown of a power plant from 
full power), then linearization must be carried out at several operating points and the 
closed loop control could be made piecewise linear by adopting the concept of gain 
scheduling. Should this concept prove to be inadequate, more advanced techniques of 
(model reference or self tuning) adaptive control [Goodwin and Sin (1984)] and 
reconfigurable control [Stengel (1991)] need to be considered. 

A block diagram of the closed loop damage control and decision system is 
proposed in Figure A.1 where [uff], [xff], {v^} and { v^} represent the sequences of 
the plant input, plant state, accumulated damage, and damage rate, respectively, 
generated as an optimal solution of the open-loop control problem via nonlinear 
programming, and [yff] is the resulting sequence of plant output, which is obtained as a 
nonlinear function g(.) of the plant state vector xff. As mentioned above, a feedback 
controller is necessary to compensate for the plant modeling uncertainties and 
disturbances such that die trajectory of the actual plant output y should be close to that 
of the desired plant output y“ which serves as the reference trajectory. The resulting 
error in the plant output is an input to the state estimator, and the feedback control 
signal ufb compensates for errors resulting from plant disturbances and sensor noise. 
The estimated state x which is obtained as the sum of the estimated state error and the 
reference state xff is fed to the structural model. The output of the structural model is 
the load vector q which contains stress, strain and other information necessary for 
damage assessment. Some of the elements of the load vector q (e.g., strain and 
temperature at the critical points) may be directly measurable as indicated in Figure A.1 
by additional measurements y. 

The closed-loop control system is partitioned into two modules. The state 
feedback control law in the first module could be formulated by using the established 
techniques of robust multi-input multi-output (MIMO) control synthesis, which rely on 
approximation of the plant dynamics by a linear time-invariant model (e.g., H 2 -based 
LQG/LTR [Stein and Athans (1987)], H 2 /H 00 optimization [Doyle et al. (1989)] or (i- 
synthesis [Doyle (1982); Doyle et al. (1982)]). However, if the plant dynamics cannot 
be approximated by piecewise linearization or time-averaging of the varying 
parameters, then selection of the control synthesis technique will depend upon the 
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specific application. The second module is a nonlinear controller which contributes to 
damage reduction in the critical plant components under feedback control and also 
serves to generate early warnings and prognoses of impending failures of the critical 
plant components. Since the damage rate v in the closed-loop control system may 
violate the specified constraints due to the additional compensation generated by the 
feedback control effort u^, a nonlinear controller is incorporated into the system to 
reduce the damage rate v. This is accomplished by corrections x cor in the reference 
trajectory and ucor in the control effort The rationale for modifying the reference 
trajectory {xff} is that, due to plant disturbances and sensor noise, it may not be 
possible to follow this trajectory without violating the damage constraints. The 
additional control effort ucor i s intended to provide a fast corrective action to the 
control input uk whenever necessary. The nonlinear control system in the outer loop, as 
shown in Figure A.1, serves two purposes, namely, (i) trimming of the linear feedback 
control signal to maintain the damage rate or accumulated damage rate vector within 
the limits, and (ii) modification of the tracking signal (i.e., the reference signal) to 
circumvent the problem of exceeding the damage rate limits, which may result from 
plant modeling errors, uncertainties and disturbances. A possible approach to synthesis 
of the nonlinear damage controller is first to postulate a mathematical structure of the 
controller and then optimize the controller parameters relative to a cost functional that 
would penalize: 

• Plant state and damage rate over the task period; 

• Final plant state and the accumulated damage. 

The outputs of the damage controller, namely x cor and u cor in Figure A. 1, need to 
be constrained to be norm-bounded to assure the system stability; xcor and ucor can 
also be considered as exogenous inputs to the linear robust control system in the inner 
loop. Therefore, the H>» bounds on xcor and u cor can be fine-tuned to satisfy the 
specified requirements of performance and stability robustness in the p-synthesis 
procedure. 


Since the damage model is highly nonlinear and is very sensitive to 
modeling uncertainties [Ray et al. (1993a); (1993b)], a combination of adaptive control 
and linear robust control techniques [Kidd (1991)] is also a viable option. However, 
such techniques are often restricted to single-input single-output (SISO) systems. The 
damage mitigation problem is MIMO because of the multiple phenomena (e.g., fatigue, 
corrosion, and creep) that may simultaneously occur at several critical points of the 
structure. Therefore, both systems-theoretic and heuristic techniques should be 
explored for synthesizing the nonlinear controller. 
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Figure A.1 Closed-Loop Damage Control System 








APPENDIX B 

RAINFLOW CYCLE COUNTING METHOD 


The load history in a real application is usually of varying amplitude and 
frequency. The fatigue damage analysis due to spectral loading is more complicated 
than that due to the constant amplitude of load. The purpose of cycle counting is to 
identify the cycles out of the spectral load history so that the strain-life data obtained 
from the experiments can be applied to the fatigue life prediction. The cycle counting 
methods, such as range-pair, rainflow, and racetrack are available in the literature 
[Fuchs and Stephens (1980)]. Among those methods, rainflow has been shown to be 
superior and yields the best fatigue-life prediction [Dowling (1983)]. The rainflow 
counting method which has been adopted in this report is briefly described below. 

Figure B.l shows the application of rainflow counting on a spectral load history. 
The stress-time curve is plotted with stress as the abscissa and time as the ordinate. 
Imagine that a raindrop is released at each peak point, flows along the direction of 
stress reversal, and then drops down. Starting from the first peak point and injecting 
the raindrops in a timed sequence, the following rules can be applied: 

(1) If the raindrop is released from the left peak point and flows to the right, the flow 
stops when it comes to a point where there is another peak point in the same 
horizontal level which is farther left than, or at the same stress level with the 
starting point 

(2) If the raindrop is released from the right peak point and flows to the left, the flow 
stops when it comes to a point where there is another peak point in the same 
horizontal level which is farther right than, or at the same stress level with the 
starting point 

(3) The flow stops when it is blocked by a previous flow path. 

In Figure B.l, the first rule is applied for the flows starting at point 2, 4 and 6. For 
example, the flow 2-3-4' stops at the point 4' since 4' is on the left of the starting point 
2. Similarly, the second rule is applied for the flows starting from the points 1 and 5. 
The flows from 3 and 7 stop at 2' and 6' due to the third rule. As can be seen from 
Figure B.l, a small cycle, 6-7-6', is extracted from a bigger cycle, 5-8-5’. After a cycle 
is closed, it can be discarded from the load history. Similariy, the cycles, 5-8-5' and 2- 
3-2', are extracted from the biggest cycle, 1-4-1'. Therefore, the large cycles in the load 
history are always retained and those cycles are most significant for fatigue damage 
prediction. The small cycles could be treated as the intermediate interruption. 





APPENDIX C 


MODELING OF PLANT DYNAMICS 

The plant under consideration in this report has been represented by a simplified 
model of a bipropellant rocket engine as shown in Figure 5.1. The physical process of 
the rocket engine consists of distributed dynamic elements which are approximately 
represented by nonlinear ordinary differential equations with time as the independent 
variable. A model solution diagram [Ray et al. (1980)] for the bipropellant rocket 
engine is shown in Figure C.l. Each block represents a physical plant component or a 
group of components. The lines interconnecting the blocks indicate the direction of 
information flow. The inputs to the plant are the commands to two oxygen valves, 
which control die liquid oxygen flowing into the prebumer and the main combustor. 



x: Torque 
ft: Angular speed 

Figure C. 1 Model Solution Diagram of a Bipropellant Rocket Engine 
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The two controlled outputs in Figure C.l are the main combustor pressure and 
the ratio of the oxygen flow and hydrogen flow. Other outputs such as turbine shaft 
speed and torque are also required to predict the damage occurring at the turbine blade. 
This model consists of eight state variables, namely turbine shaft speed, pump (LH 2 - 
fuel) mass flow rate, prebumer gas pressure, prebumer gas density, combustor gas 
pressure, combustor gas density, and the two flow rates of oxidant 

The model equations are formulated by using approximation of the integral 
forms of the conservation relations for mass, momentum, and energy, semi-empirical 
relationships for fluid flow and heat transfer, and state relations for thermodynamic 
properties of the working fluids. Model parameters are evaluated from the design and 
experimental data of the manufacturer [NASA (1991)]. Major assumptions in addition 
to the lumped parameter approximation of the model equations are: 

(1) Uniform fluid flow over the pipe cross-section; 

(2) No heat gain or loss for the fluid flowing through the pipes; 

(3) Perfect thermal insulation between plant components and the environment; 

(4) Negligible pressure drop due to velocity and gravitational heads in the fluid 
paths; 

(5) Constant density of the oxygen and hydrogen in the liquid phase; 

(6) Validity of perfect gas law inside the prebumer and main combustor; 

(7) Choked flow through the turbine and nozzle; and 

(8) Representation of the valve dynamics by a first order lag. 

The following subsections formulate the governing equations of the plant 
components including the liquid hydrogen pump, pump turbine, liquid hydrogen 
pipeline, two oxygen inlet valves, prebumer and main combustor as shown in Figure 
Cl. 


The control input, measured output, and plant state variables are listed in Table 
C.1 along with their steady-state values at two operating conditions. The system 
matrices and the eigenvalues of the linearized state space model at two steady-state 
operating conditions are listed in Table C.2. 

Pump Model: 

The pump driven by the turbine feeds pressurized liquid hydrogen into the 
prebumer. The empirical relation describing the pump head-flow-speed characteristics 
[Hicks and Edwards (1971); Ray (1976)] is given as 


AP 


prop _ 


Ph^pmp^ 


W 


prop 


V 


Ph«, 


propy 


(C.1) 


where AP pmp is the pressure increase between the inlet and outlet of the pump, W pmp 
and Ph are the mass flow rate and mass density of liquid hydrogen, respectively, and 
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Qpmp is the pump speed. Eq. (C.l) can be normalized with respect to a normal 
operating point (which is characterized by the values of the process variables marked by 
the superscript *) such that 


APpmp ! (Ph^pmp^) 


AP*/(pV 2 ) 


= f 


Wpmp / Ph^pmp 

5 * t 

l W /p Q 


(C.2) 


Assuming that the density of liquid hydrogen remains constant within the operating 
range, the pump characteristics described in eq. (C.2) can be approximated by a second 
order polynomial based on the manufacturer's data [NASA (1991)] as: 
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The mass flow rate, W p pp, of liquid hydrogen and the angular speed, Q pm p, of the 
turbo-pump are state variables. The other coefficients, ao, ai, a 2 and the variables at 
normal operating point are given as pump data. Then, the pressure head, APp mp , and 
the pump outlet pressure, Pdei, are computed by 

APpmp = AP*(a pmp /Q*)2y (C.4) 

Pdel = APpmp + Psuc (C.5) 

where the pump inlet pressure, Psuc, is assumed to be constant and equal to the pressure 
at the source of liquid hydrogen. Once APpmp is known, the output power of the pump 
can be obtained as: 


Output power = AP pmp - JS E (C.6) 

The input power and the torque T pm p required to drive the pump are then computed 
from the pump efficiency Tipmp as: 

Input power = AP pmp J* and 'tpmp — AP pm p — — SEE (C.7) 

"h T ‘pmp "h^lpmp^^pmp 

where the efficiency is also approximated by a second order polynomial based on the 
manufacturer’s data [NASA (1991)] as: 
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(C.8) 


Assuming perfect insulation between the pump and the environment, the energy loss is 
completely transformed into heat and causes the increase of enthalpy (and temperature) 
for the liquid hydrogen as presented below: 
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where J is the conversion factor from the thermal unit to the work unit, and C p h is the 
specific heat of liquid hydrogen at constant pressure. 


Pipeline Model for Liquid Hydrogen Flow: 


For the liquid hydrogen flowing over the pipeline, the head loss occurs due to 
friction along the wall of the pipeline. Let A and / be the average area of cross section 
and the length of the pipeline respectively. Then, the pipeline model is formulated by 
approximating the integral form of the momentum equation as: 


^ (Ph'Av ) — P dg]A P pb r A F 
^ ~“(Ph Av ) = f*del ” Ppbr “ APfjjc 

=» {^“Pdd-Ppta-APftic < Clt » 

where P p br is the pressure of the prebumer and the frictional pressure drop is obtained 
as [Shames (1962)]: 


I WpjniJ W™,,, 

APfric = Kfrie P°*P PSE (C.11) 

Ph 

The loss of heat from the pipe line is assumed to be balanced by the gain of heat due to 
friction. Thus, the prebumer inlet temperature is taken to be equal to the pump outlet 
temperature. 

Oxygen Inlet Valve Models: 

Two oxygen valves control the oxidant flow into the prebumer and combustor. 
Let P 0 be the pressure inside the oxygen tank. The mass flow rates of oxygen are 
computed from the orifice equation [Blackburn et aL (I960)] as 

Wpfor — Kpjjf £pbr-\/(Po ~ Ppbr )Po (C.12) 

Wcmb “ Kcmb ^cmbV^o ~ Pcmb)Po (C.13) 

where p 0 is the density of liquid oxygen, and ^ p br and £c m b are the actual valve 
positions for the prebumer and combustor respectively. The dynamics of the valves 
subject to the command inputs are assumed to behave as the first order lag such that 
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d^pbr _ £pbr u pbr 
dt ®pbr 


(C.14) 


(C.15) 


d^cmb _ £cmb u cmb 
dt ®cmb 

where u p b r and Uc m b are the electrical commands to the oxygen valves of prebumer and 
combustor, and 0pbr and Qcmb arc the time constants of the respective values. 


Turbine Model: 

For simplicity of modeling, choked flow [Blackburn et al. (I960)] is assumed 
through the turbine, and no loss of pressure and enthalpy occurs between the prebumer 
and the turbine. Therefore, die turbine inlet pressure and temperature are assumed to be 
identical to the pressure and temperature in the prebumer. Then, the mass flow rate of 
hot gas flowing through the turbine can be expressed by 

(C.16) 

V^pbr 

where Tjjbr is the absolute gas temperature in the prebumer. For an isentropic process, 
the turbine exhaust pressure, Pirb> and ideal outlet temperature, T^j, are related to the 
turbine inlet pressure, Ppbr, and inlet temperature, Tpbr, by 


Ppbr 1 ^T p br^=Pt lb 1 ^T tlb ? (C.17) 

where Y = Cp/C v for the gas. Following the information flow in the model solution 
diagram in Figure C.l, die turbine exhaust pressure, Ptrb» is approximated by the 
combustor pressure plus the head loss due to friction through the injector tubes. A 
relationship similar to eq. (C.12) or (C.13) is used to obtain the mass flow and pressure 
drop in the injector (which connects the turbine exhaust with the combustor chamber): 

Wtri, = ^ubf V^trb ~ Pcmb )Pcmb (C.l 8) 

The above eq. (C.l 8) has been approximated by using the density p C mb of the 
combustor gas instead of the average gas density in the injector. The rationale for this 
approximation is that eqs. (C.16), (C.17) and (C.18) can be combined to obtain the 
turbine exhaust pressure Ptrbe via a closed form relation, and the ideal oudet 
temperature can be direcdy evaluated from eq. (C.17). The actual oudet temperature is 
estimated from the turbine efficiency, which is given by the following empirical 
equation [NASA (1991); Ray (1976)]: 
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where tj*, Q* and AH* are the turbine efficiency, speed and ideal enthalpy drop at the 
designed operating condition. The turbine speed, Qtrt» is identical to the pump speed, 
and therefore is a state variable. The ideal enthalpy drop, AH^i, is given by 


AHtrbi — Cppbr (Tpbr “ Tpbi) (C.20) 

The actual enthalpy drop is then obtained as AHtrb = AHubi Tltrb- Since the difference 
between the ideal enthalpy drop and the actual enthalpy drop is converted to irreversible 
loss of energy and no other loss is assumed to be prevalent, the turbine torque Ttrb is 
obtained as 

%b — WtibAHtibA^pmp (C.21) 

where Qpmp is the angular speed of the turbo-pump. Having known ipb, the equation 
of motion for the turbine is given from conservation of angular momentum [Shames 
(1962)]: 


A pmp “ T trb t pmp 

where I pm p is the moment of inertia of the turbo-pump assembly. 
Combustor Model: 


(C.22) 


The combustor model is formulated via approximation of the integral form of 
conservation of mass flow and internal energy. Assuming the gas inside the main 
combustor chamber to be perfect, the variables describing the thermodynamic states in 
the combustor chamber are related by the ideal gas law [Zemansky and Van Ness 
(1966)]: 


P cmb - R Pcmb Tcmb (C.23) 

where R is the characteristic gas constant In the model, the pressure, P C mb> and the 
density, Pcmb. are selected as the state variables. The absolute temperature, Tcmb> of the 
combustor gas is computed from eq. (C.23) as a function of these two state variables. 
The exhaust hot gas, which generates the thrust for the rocket engine, is obtained as the 
choked flow through the nozzle. Then, the mass flow leaving the combustor can be 
computed as 
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Having known the nozzle exhaust flow from eq. (C.24), the gas density inside the 
combustor chamber, Pcmb> is computed from the conservation of mass such that 

dPcmb _ Wtib + W cm b — W noz 
Vcmb 


dt 


( 025 ) 
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where V cm b is the volume of the main combustor chamber. Assuming prefect 
insulation without any energy loss to the environment, the conservation of internal 
energy yields: 

— (C v VPT) cmb = WjjbH^b + WgnjbCpQTQ — W noz Cp Cmb T cin | ) + (C.26) 

where the last term, CfueiW cm b, is the heat generated from consumption of the oxidant 
in the fuel-rich environment Substituting pT by P/R from the ideal gas law, eq. (C.26) 
can be written as 


dPqnb _ 
dt 
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R 


(C.27) 


Prebumer Model: 


The prebumer model is formulated similar to the combustor model with the gas 
density, p p b r , and the pressure, Pnbr* as the state variables. The absolute gas 
temperature is computed from the ideal gas law. From the conservation of mass and 
energy, two equations similar to eqs. (C.25) and (C.27) are given as 
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Table C.1 Plant Control Input, Measured Output and State Variables 

Steady-state values for combustor pressure 


Input variables: 

2700 psi 

3000 psi 

ui : Command to the prebumer O 2 valve 
(dimensionless) 

0.655466 

0.838256 

U 2 : Command to the combustor O 2 valve 
(dimensionless) 

0.697910 

0.807391 

Measured output variables: 

yi : Combustor gas pressure (psi) 

2700 

3000 

y 2 : 02 /H 2 ratio 

6.02 

6.02 

State variables: 

xi : Turbine (pump) shaft speed (rad/sec) 

3606.40 

3904.71 

X 2 : LH 2 -fuel mass flow rate (Ibm/sec) 

115.363 

128.206 

X 3 : Prebumer gas pressure (psi) 

3737.11 

4240.40 

X 4 : Prebumer gas density (lbm/in 3 ) 

0.3207xl0“ 3 

0.3571x10-3 

X 5 : Combustor gas pressure (psi) 

2700.00 

3000.00 

X 6 : Combustor gas density (lbm/in 3 ) 

0.1038xl0- 3 

0.1153x10-3 

X 7 : Position of the prebumer O 2 valve 
(dimensionless) 

0.655466 

0.838256 

X 8 : Position of the combustor O 2 valve 
(dimensionless) 

0.697910 

0.807391 
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Table C.2 System Matrices for the Linearized Plant Model 

x = Ax + Bu 
System equations : y = Cx + Du 


System matrices for combustor pressure at 3000 psi : 


A = 


-.1685D+02 
0.6064D+02 
0.1215D+02 
0.0000D+00 
02268D+01 
0.0000 D+00 
0.0000 D+00 
0.0000 D+00 


-.6154D+02 
-6002D+03 
0.9678D+03 
0.8333D-03 
0.0000 D+00 
0.0000 D+00 
0.0000 D+00 
0.0000 D+00 


0.4105D+02 
-2000 D+ 02 
-.1655D+04 
-3547D-04 
02210D+03 
0.5349D-05 
0.0000 D+00 
0.0000 D+00 


-.1240 D+ 09 
0.0000 D+00 
0.3953D+10 
-2382D+03 
-.9147D+09 
0.635 1D+02 
0.0000 D+00 
0.0000 D+00 


-.4582D+02 
0.0000 D+00 
0.0000 D+00 
0.0000 D+00 
-4586D+04 
-.5675D-04 
O.OOOOD+OO 
O.OOOOD+OO 


02544D+09 

0.0000D+00 

O.OOOOD+OO 

O.OOOOD+OO 

0.3122D+11 

-.8664D+03 

O.OOOOD+OO 

O.OOOOD+OO 


O.OOOOD+OO O.OOOOD+OO 
O.OOOOD+OO O.OOOOD+OO 
0.3206D+07 O.OOOOD+OO 
0.7564D-01 O.OOOOD+OO 
O.OOOOD+OO 0.8116D+07 
O.OOOOD+OO 0.1915D+00 
-.3333D+03 O.OOOOD+OO 
O.OOOOD+OO -3333D+03 


O.OOOOD+OO O.OOOOD+OO 
O.OOOOD+OO O.OOOOD+OO 
O.OOOOD+OO O.OOOOD+OO 
O.OOOOD+OO O.OOOOD+OO 
O.OOOOD+OO O.OOOOD+OO 
O.OOOOD+OO O.OOOOD+OO 
03333D+03 O.OOOOD+OO 
O.OOOOD+OO 0.3333D+03 


c = 


O.OOOOD+OO O.OOOOD+OO O.OOOOD+OO 
O.OOOOD+OO -.4682D-01 -.1443D-03 


O.OOOOD+OO 1.0000D+00 O.OOOOD+OO 
O.OOOOD+OO -.8228D-03 O.OOOOD+OO 


O.OOOOD+OO 

O.OOOOD+OO ' 

1" O.OOOOD+OO 

O.OOOOD+OO ' 

0.7080D+00 

0.6721D+01 _ 

~ L O.OOOOD+OO 

O.OOOOD+OO 


Eigenvalues of matrix A (sec* 1 ): 

Real part Imaginary part 


Real part Imaginary part 


-0.402555D+04 
-0.158369D+04 
-0.142605D+04 
-0.451 191D+03 


O.OOOOOOD+OO 

O.OOOOOOD+OO 

O.OOOOOOD+OO 

0202507D+03 


-0.451191D+03 

-0254363D+02 

-0333333D+03 

-0333333D+03 


-0.202507D+03 

O.OOOOOOD+OO 

O.OOOOOOD+OO 

O.OOOOOOD+OO 
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System matrices for combustor pressure at 2700 psi : 


A = 


-1532D+02 
0.5594D+Q2 
0.1115D+02 
O.OOOOD+OO 
0.1780D+01 
0.0000 D+00 
0.0000 D+00 
0.0000 D+00 


-.5684D+02 0.4304 D+ 02 
-.5382D+03 -2000D+02 
0.9505D+03 — .1450D+04 
0.8333D-03 -3105D-04 
0.0000D+00 0.2200 D+ 03 
O.OOOOD+OO 03399D-05 
O.OOOOD+OO 0.0000D+00 
O.OOOOD+OO O.OOOOD+OO 


-.1234 D+09 -.4825D+02 
O.OOOOD+OO O.OOOOD+OO 
0.3844D+10 O.OOOOD+OO 
-3360D+03 O.OOOOD+OO 
-.8969D+09 — .4417D+04 
0.6292D+02 -3270D-04 
O.OOOOD+OO O.OOOOD+OO 
O.OOOOD+OO O.OOOOD+OO 


O.OOOOD+OO O.OOOOD+OO 
O.OOOOD+OO O.OOOOD+OO 
03576D+07 O.OOOOD+OO 
0.8438D-01 O.OOOOD+OO 
O.OOOOD+OO O.8477D+07 
O.OOOOD+OO 0.2000 D+00 
-3333D+03 O.OOOOD+OO 
O.OOOOD+OO -.3333D+03 


B = 


O.OOOOD+OO 

O.OOOOD+OO 

O.OOOOD+OO 

O.OOOOD+OO 

O.OOOOD+OO 

O.OOOOD+OO 

03333D+03 

O.OOOOD+OO 


c = 


O.OOOOD+OO O.OOOOD+OO 0.0000D+00 O.OOOOD+OO 1.0000D+00 
O.OOOOD+OO -.5203D-01 -.1124D-03 O.OOOOD+OO -.7566D-03 


O.OOOOD+OO O.OOOOD+OO ' 
0.8777D+00 0.7801D+01 _ 


D = 


O.OOOOD+OO 

O.OOOOD+OO 


Eigenvalues of matrix A (sec* 1 ): 


Real part 

-0.386847D+04 

-0.141679D+04 

-0.139069D+04 

-0.413837D+03 


Imaginary part 

O.OOOOOOD+OO 

O.OOOOOOD+OO 

O.OOOOOOD+OO 

0340277D+03 


Real part 

-0.413837D+03 

-0.1%700D+02 

-0333333D+G3 

-0333333D+G3 


0.2618D+09 

O.OOOOD+OO 

O.OOOOD+OO 

O.OOOOD+OO 

03125D+11 

-.8665D+03 

O.OOOOD+OO 

O.OOOOD+OO 

O.OOOOD+OO ‘ 

O.OOOOD+OO 

O.OOOOD+OO 

O.OOOOD+OO 

O.OOOOD+OO 

O.OOOOD+OO 

O.OOOOD+OO 

03333D+03 . 

O.OOOOD+OO 

O.OOOOD+OO 

O.OOOOD+OO 1 
O.OOOOD+OO J 


Imaginary part 

-0.240277D+03 

O.OOOOOOD+OO 

O.OOOOOOD+OO 

O.OOOOOOD+OO 



APPENDIX D 


MODELING OF STRUCTURAL DYNAMICS 


One of the critical points where severe damage is likely to occur is assumed to 
be located at the root of die turbine blades. Since both the damage and subsequent 
crack initiation sites are confined within a small region of the blade, a linear elastic 
approach is adopted for macroscopic modeling of the structural dynamics and to predict 
transient stresses at the point of potential failure. In this approach, the blade geometry, 
properties of the blade material, and state variables of the plant dynamic model (see 
Appendix C) are used as inputs to the finite element analysis program to generate a 
discretized representation of the blade structure and its loading conditions. The 
resulting stiffness matrix, mass matrix, and force vector are then used to obtain a modal 
solution for the displacements. In the last step, the stress-displacement relations from 
the finite-element analysis are used to predict the stresses at the critical point(s) of the 
blade structure. The general approach is shown in the flow chart of Figure D.l. 
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Figure D. 1 Flow Chart for the Structural Model of a Turbine Blade 
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For simplicity of computation, the turbine blade is approximated to behave like 
a beam as shown in Figure D.2. Let the torque applied to the turbine shaft at time t be 

T(t). Then, the (time-dependent) mean load, P m (t), per unit length exerted on a typical 
blade is computed by 




T(t) 

NfiLim 


(D.l) 


where Nb is the number of blades of the turbine, L is the length of the turbine blade and 
r m is the mean radius of the blade for rotation about the shaft axis. The mean load, 
P m (t), serves as the static load on the blade at time L The dynamic load is assumed to 
oscillate with an amplitude proportional to the static load and a frequency equal to the 
product of the shaft speed, Q(t), and the number, Ny, of stationary vanes. The total 
load per unit length acting on one turbine blade at time t is then given as the sum of the 
static and dynamic loads per unit length as: 

P(t) = Pm(t) + PA(t) sm(Q(t)N v t) (D.2) 

where Pa( 0 = P P m (t) and the proportionality constant p is assumed to be an 
approximate representation of steady-state fluid dynamics within the turbine. 


Figure D.2 Loading of a Turbine Blade 



With the total load per unit length given by eq. (D.2), the finite element method 
is used to estimate the stress and strain at the root of the turbine blade. The blade is 
modeled to consist of two beam elements with three nodes equally spaced as shown in 
Figure D.3. The model has a total of eighteen degrees of freedom where each node has 
six degrees of freedom with three translational and three rotational modes. The first 
node is fixed and thereby all six degrees of freedom associated with the first node are 
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suppressed. This is achieved via the penalty approach [Chandrupatla (1991)] in which 
the spring stiffness at each of the six degrees of freedom at the first node is set to a very 
large value. In contrast, the use of the elimination approach [Chandrupatla (1991)] 
would reduce the order of the model to twelve. Ease of implementation in the existing 
finite-element analysis code is the rationale for using the penalty approach as opposed 
to the elimination approach. 



Figure D.3 Finite Element Model of a Turbine Blade 

Let a represent the nodal displacements of a beam element, where a is a 12x1 
vector with 6 degrees of freedom at each node, the displacement field u(x) along the 
element can be expressed in terms of the nodal displacements a and an interpolation 
matrix N: 


u = N a (D.3) 

where u(x) contains three translational and three rotational displacements at location x 
of the beam element, and each element of the matrix N is a polynomial interpolation 
function of x. The dimension of u is 6x1, and the dimension of N is 6x12. The 
generalized strain field e is then obtained via differentiation of the displacement field 
as: 


e = Lu = LNa = Ba (D.4) 

where L is a differential operator and B = L N. The dimensions of L and B are 4x6 and 
4x12, respectively. The generalized strain vector e contains 4 elements representing the 
generalized strains due to twisting, stretching and bending in two directions orthogonal 
to the beam axis. The twisting action produces shear stresses. The stretching and 
bending actions produce tensile or compressive stresses. For linear elastic materials, 
the generalized stress vector o is related to the generalized strain vector e by a matrix 
D such that 


a = D e (D.5) 

The next step is to model the beam dynamics with nodal displacements as the 
variables of motion. From eq. (D.3), the total kinetic energy, T, is computed by 
integrating the kinetic energy of the mass, pdx, of the differential element along the 
entire length, / , of the element : 
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T = ij' fiTpudx 

= \\‘ a a T N T pNadx 

= i iT (/o N T P Ndx )a 

= — a T Ma 
2 


(D.6) 


where p is the diagonal matrix containing the mass density terms and M = 
N T pN dx is generally called the mass matrix of the element The potential energy, 
V, can also be computed by integrating the strain energy of the element as given below: 


v = Ij'e T odx 

= ij' e T De dx 

= -f! a T B T DBa dx 
2 JO 

= r T (lo B T DBdx)a 

= — a T Ka CD-7) 

2 


where K = J* B T DB dx is the stiffness matrix. Since the energy is additive, the mass 

and stiffness matrices of the whole structure consisting of the finite element meshes can 
be obtained by the principle of superposition. The matrices of the individual elements 
are added together at each degree of freedom in the formulation of matrices for the 
entire structure. The matrices M and K are both real symmetric since p and D are 
symmetric. Furthermore, die mass matrix M is positive definite because the kinetic 
energy is always positive with a nonzero velocity, a . The stiffness matrix K is positive 
semidefinite due to nonnegative strain energy. From the Lagrange equations of motion, 
the beam dynamics in the absence of any damping is derived below : 

d f 3(T- V) ^ 3(T- V) _ f 
dt^ 3a J 3a 

=> Ma + Ka = f (D.8) 


where f is the vector of the concentrated nodal forces applied to the degrees of freedom. 
If the distributed load, p, is applied, the work done by p along the length of the beam 
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should be equal to the work done by the equivalent concentrated forces, f, on the nodal 
displacements such that 

a T f = Jq u T p dx = Jq a T N T p dx = a T J^ N T p dx (D.9) 

Since eq. (D.9) is true for all nodal displacements a, the equivalent forces can be 
calculated as: 


f = Jq N T p dx (D.10) 

Eq. (D.8) does not include the forces caused by damping. If the damping effects 
need to be considered, the damping coefficients can be approximated as a linear 
combination of the mass and stiffness matrices [Weaver (1991)] such that C = aM + 
pK. Then, eq. (D.8) is modified as 

Ma + (aM+pK)a + Ka = f (D.ll) 

Since eq. (D.ll) contains a set of coupled differential equations, its numerical solution 
through direct integration is time-consuming. An alternative approach is to find a 
transformation matrix and diagonalize die coupled equations into a set of decoupled 
equations that have the same form as ordinary second order differential equations. This 
set of decoupled equations can be easily solved and transformed back to the original 
nodal displacements. If the nodal forces, f, are in the form of general periodic 
functions, then they can also be expanded in Fourier series containing trigonometric 
functions such that the exact solution can be derived. Another advantage of this 
approach is that some of the decoupled equations associated with high natural 
frequencies lead to small amplitudes of vibrations and may therefore be neglected in 
computation. Consequently, computation time is saved without any significant loss of 
accuracy. 

To decouple eq. (D.ll), first consider the following eigenvalue problem: 

Kx = © 2 Mx (D.12) 

Since eqs. (D.6) and (D.7) yield the symmetric positive definite mass matrix M and the 
symmetric positive semidefinite stiffness matrix K, the eigenvalues associated with eq. 
(D.12) are all real and nonnegative. Let ©i 2 , CO 2 2 , ••• , ©n 2 be the eigenvalues of eq. 
(D.12) and xi, X 2 , x n be the corresponding eigenvectors, where ©j and xj are 
generally referred to the natural frequency and mode shape associated with the i* mode 
of the structure [Weaver (1991)]. For any two distinct eigenvalues, ©i 2 and ©j 2 , eq. 
(D.12) can be written as 


Kxi = ©i 2 Mxi 
Kxj = ©j 2 Mxj 


(D.13a) 

(D.13b) 
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Premuldplying both sides of eq. (D.13a) by Xj T and eq. (D.13b) by Xj T , it becomes: 

XjTKxj = fi)i 2 Xj T Mxi (D.14a) 

XjTKxj = ©j 2 xjTM Xj (D.14b) 

Since M and K are both symmetric, Xj T K Xi = Xi T K Xj and xj T M Xi = xi T M Xj. 
Subtracting (D.14b) from (D.14a) results m 


0 = ( ©i 2 - ©j 2 ) xjTM xi (D.15) 

If ©i 2 * ©j 2 , then both XjTMxi and XjTKxi must be equal to zero. This result leads to 
the following important property: if all the eigenvalues are distinct, the eigenvectors are 
mutually orthogonal with respect to the mass and stiffness matrices. In fact, if an 
eigenvalue Xj for real symmetric matrices M and K has multiplicity m, then there exist 
exactly m linearly independent eigenvectors corresponding to these m eigenvalues 
[Meirovitch (1980)]. It is always possible to choose such combinations that these 
eigenvectors are mutually orthogonal and orthogonal to eigenvectors belonging to the 
remaining eigenvalues. Let each eigenvector, xi, be normalized with respect to the 
mass matrix M, i.e., XjTMxi =1 for i = 1,2, ••• n. Then, following either eq. (D.14a) or 
(D.14b), it can be deduced that XiTKxj = ©i 2 . The above results can be written as: 


Xj T M Xj = 0 ifi.*j Xj T M Xj = 1 if i = j (D.16a) 

Xi T Kxj =0 if i*j Xi T Kxj = ©i 2 if i=j (D.16b) 


Let 'P = [ xi X2 ••• x n ]and A = 


©!* 0 
0 © 2 2 


0 


0 


© 


n J 


Then, eqs (D.16a) and 


(D.16b) can be expressed as 

^¥ = 1 and ^ T K 1 P=A (D.17) 

Then, with premultiplication by *PT, eq. (D.ll) can be decoupled as: 


q + (ctI+PA)q + Aq = Y T f (D.18) 

where q=Y _1 a and exists because of linear independence of eigenvectors, 
Xi, i = 1,2, ••• n. After solving eq. (D.18) for the vector q, the nodal displacement 
vector can be computed from the transformation a = Yq. The generalized strains and 
stresses at the selected points are then recovered from eqs. (D.4) and (D.5). 
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Due to the damping effects, the individual components of the displacement 
vector, a, may have different phase angles. It is expensive to compute the stress at a 
certain point by simple superposition of the effects of twisting, stretching and bending. 
It is assumed in the simulation that the damping effects are negligible and all 
components of the displacement vector have an identical phase angle. The stresses 
caused by stretching and bending are added together for the normal stress. The shear 
stress is obtained from twisting. The maximum principal stress is an approximate 
representation of the stress at the comers of the root of a turbine blade, which is used to 
estimate the fatigue strain-induced damage rate. 



APPENDIX E 


MODELING OF INTEGRATED FATIGUE-CORROSION-CREEP DAMAGE 

Experimental results indicate that, in gas turbine materials such as Ni-base 
superalloys, reduction of fatigue life at elevated temperatures (as low as 30% of the 
melting point) is significantly influenced by the gaseous environment [Courtney (1990); 
Suresh (1991)]. The contributions of creep and corrosion on structural damage become 
significant at elevated temperatures due to the effects on micromechanisms of fatigue 
fracture. For example, the fatigue failure at a low temperature usually results from 
transgranular fracture but it may occur due to intergranular fracture at an elevated 
temperature. The impact of elevated temperatures on damage could be both beneficial 
and detrimental. Beneficial effects include a dispersal of the slip, and the detrimental 
effects are usually caused by creep and corrosion. In fact, the fatigue life at a high 
temperature may increase with an increase in the frequency of strain reversals [Suresh 
(1991)]. While the fatigue damage is cycle-dependent, the creep damage and corrosion 
damage are time-dependent 

Although the time-dependent part of the fatigue-creep damage is thermally 
activated and results from both creep and corrosion, linear elastic fracture mechanics 
can be used for damage prediction provided that the stress intensity factor range 
remains the controlling parameter for crack propagation. This is valid if the zone of 
inelastic deformation at the crack tip is small in size compared to the uncracked length. 
However, at very high temperatures and low cycles, the elastic fracture mechanism may 
not be prevalent as the inelastic zone becomes appreciably large. Furthermore, the 
material properties such as the tensile strength and fracture toughness change 
significantly at elevated temperatures relative to their nominal values, Tlie dynamics of 
creep damage may be modeled from these perspectives. The total strain rate e 
consisting of elastic strain rate e e and plastic strain rate e p is modeled for an elastic- 

nonlinear viscous material under uniaxial tension in terms of the instantaneous stress o 
and the stress rate a as [Lemaitre (1992)]: 


e = e c +e p 


where 



and 



\ n c 


v a yj 


(E.1) 


where the material parameters are: Young's modulus E; yield strength Gy; reference 
creep strain rate e y ; and power law creep exponent n c . Eq. (6) may have to be 

modified for multi-axial stress conditions in the critical components. In that case, the 
effective stress needs to be defined as a function of the stress tensor. 

Some investigators, for example [Saxena (1988)] and [Zamrik and Davis 
(1990)], have computed the total damage as the sum of corrosion-fatigue and creep 
damages that are obtained independently. That is, the cumulative damage Dcf in creep- 
fatigue cycling has been defined upon completion of n cycles as: 


Dcf = Df + Dc 


(E.2) 
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where Dc is the creep damage part; Df is the fatigue damage part; and Dcf is 
normalized to unity at the onset of failure. Similar to the linear damage rule for 
cumulative fatigue damage, the creep damage is often based on time summation as: 




(E.3) 


where Nfh is the number of cycles to failure in creep-fatigue; tr is the time to creep- 
rupture; and th is the hold time. However, the above time summation model lacks the 
input of the material properties. Based on observed failure modes for the creep-fatigue 
interaction, ductility is considered to have a major influence on failure of engineering 
materials. Zamrik and Davis (1990) have proposed a model of the creep damage part 
D c based on the ductility exhaustion concept as follows: 


n 

th 

A£~, 

* 

D c =£ 

N*J 

cp 

^ / - * \ 

dt 

i=l 

0 

k 

*(^cp). 



where ACq, is the rate of creep strain range during the hold-time period, and < d(AE C p) 
is the material creep ductility as a function of Ae^,. If the creep strain accumulation 
during the hold-time period is assumed to obey the power law, then 

A£cp = A(Aa) q t r_1 (E.5) 

where A, r and q are constants, and Ao is the stress range at time t during hold period. 
Pineau (1989) also proposed a ductility exhaustion model where the creep damage D c 

is expressed in terms of the instantaneous values of plastic strain rate e p , stress a, and 
temperature T as: 


D c (t) = 


\( e p (a(x),T(x))^ ^ 

^e s (o(x),T(x)) J x t r (o(x),T(x)) 


(E.6) 


where t is the current time, and tr is the total time to rupture due to creep, and e s is the 
stationary creep rate [Coutsouradis et al. (1978)]. The two models in eqs. (E.5) and 
(E.6) in view of the fact that ductility is the primary parameter for creep damage at 
elevated temperatures. 

Modeling of corrosion-fatigue crack growth in gaseous environments has been 
reported by several investigators including Wei (1989). The governing equations for 
damage dynamics are formulated on the basis of mass balance and reaction kinetics 
where the state variables are gas pressure at the crack tip and percent of surface 
reactions. In general, the corrosion-fatigue crack growth rate (da/dN)cf can be obtained 
as a modification of the pure fatigue growth rate (da/dN)f by including the states of gas 
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pressure and surface reactions. If the dynamics of corrosion are likely to be slow 
relative to fatigue dynamics except for operations at elevated temperatures. Depending 
on the operating conditions, it might be possible to e limin ate the above state variables 
by varying the parameters of the fatigue crack growth equation slowly with time. 
Therefore, the constraint parameters p and r in eqs. (2.4) and (2.5) may also be chosen 
as slowly varying functions of time to incorporate the effects of corrosion. 

Following the procedures described in Section 4, the cycle-dependent corrosion- 
fatigue crack growth rate (da/dN) C f can be converted to the time-dependent crack 
growth rate (da/dt) C f. The corrosion-fatigue damage is defined, similar to eq. 4.31, as 

D C f = a/a* where a* is the critical crack length. Then, following the structure of eq. 
(2.2), a dynamic model of the combined effects of creep and fatigue damage at a single 
critical point is proposed as follows: 


dD c f 



= MDcKt), D CT (t), q(x, t), t); 

= hcr(Dcf(t), Dcr(t), q(x, t), t); 


DcfOo) =D c fo; 
Dcr(to) =DciOi 


hcf>0 Vt 
ho- > 0 Vt 


(E.7) 

(E.8) 


In the above equation set, the damage vector is expressed as v(t)=[D c f D C ] T . 
For multiple critical points, the damage vector includes the damage components at all 
points. Finally, a scalar measure D of total damage is expressed as a function of the 
corrosion-fatigue damage and creep damage, and the end of service life occurs at D=l. 
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